Finite-size scaling and boundary effects in two-dimensional valence-bond-solids 
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Various lattice geometries and boundary conditions are used to investigate valence-bond-solid 
(VBS) ordering in the ground state of an 5* = 1/2 square-lattice quantum spin model — the J-Q 
model, in which four- or six-spin interactions Q are added to the standard Heisenberg exchange 
J. Ground state results for finite systems (with up to thousands of spins) are obtained using an 
unbiased projector quantum Monte Carlo method. It is found that great care has to be taken when 
extrapolating the order parameter to infinite lattice size, in particular in cylinder geometry. Even 
though strong VBS order exists in two dimensions, and is established clearly with increasing system 
size on LxL lattices (or x Ly lattices with a fixed aspect ratio L^/Ly of order 1), only short-range 
VBS correlations are observed on long cylinders (when La; — >■ cxj at fixed Ly). The correlation length 
increases with the cylinder width, until long-range order sets in at a "critical" width. This width is 
very large even when the 2D order is relatively strong. For example, for a system in which the order 
parameter is 70% of the largest possible value, Ly = 8 is required for ordering. Extrapolations of the 
VBS order parameter based on correlation functions (the square of the order parameter) for small 
LxL lattices can also be misleading. For a 20%-ordered system results for L up to ~ 20 appear 
to extrapolate clearly to a vanishing order parameter, while for larger lattices the scaling behavior 
crosses over and extrapolates to a non-zero value (with exponentially small finite size corrections). 
The VBS order parameter also exhibits interesting edge effects related the known emergent U(l) 
symmetry close to a "deconfined" critical point, which, if not considered properly, can lead to wrong 
conclusions for the thermodynamic limit. The observed finite-size behavior for small LxL lattices 
and long cylinders is very similar to that predicted for a Z2 spin liquid. The results therefore 
raise concerns about recent numerical work claiming Z2 spin liquid ground states in 2D frustrated 
quantum spin systems, in particular, the Heisenberg model with nearest and next-nearest-neighbor 
couplings. Based on the results presented so far, a VBS state in this system cannot be ruled out. 
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PACS numbers: 75.10.Kt, 75.10.Jm, 75.40.Cx, 75.40.Mg 



I. INTRODUCTION 



A valence-bond solid (VBS) is a state of a quantum 
spin system in which there is no magnetic long-range 
order, but lattice symmetries (translational and some 
times rotational) are broken due to dimerization or, more 
generally, polymerization of the system into one with 
a larger unit cell than the underlying lattice. One can 
think of the spins within a unit cell of a VBS (or within 
different groups of spins in a large complex unit cell) 
as having an enhanced probability of forming a total 
spin singlet. In this paper, manifestations of VBS order 
in ground states of finite systems are investigated, us- 
ing unbiased quantum Monte Carlo (QMC) simulations 
of 5 = 1/2 spins on the two-dimensional (2D) square 
lattice with interactions — Heisenberg exchange supple- 
mented by certain multi-spin interactions — leading to 
columnar order in the thermodynamic limit ii The ap- 
proach to the infinite-size 2D limit is investigated for 
different boundary conditions. The models considered 
can be tuned from strong to weak VBS order (and also 
through a critical point), enabling bench-mark investiga- 
tions of asymptotics and cross-over behaviors. In partic- 
ular, consequences of near-criticality of the VBS order on 
the finite-size behavior can be examined in detail. The 
stability of VBS order on long cylinders {L^ x Ly lat- 
tices with Lx ^ Ly) is also addressed. This geometry 
is often used in density matrix renormalization (DMRG) 



studies^i^ with recent intriguing results pointing to the 
absence of VBS order and the existence of spin liquids in 
frustrated models whose ground states have been debated 
for a long time^Ti^ 

In the following introductory sections, several back- 
ground facts motivating further studies of VBS order are 
discussed and some of the known properties of VBS states 
are briefly reviewed. The purposes of the studies reported 
here will then be detailed, followed by an outline of rest 
of the paper. 



A. VBS states and frustrated interactions 

VBS states have been known for a long time to exist 
in ID frustrated quantum spin chains. In particular, in 
the 5 = 1/2 Heisenberg chain with nearest- and next- 
nearest- neighbor couplings Ji and J2, the ground state 
at coupling ratio g = J2/ Ji = 1/2 is exactly a product 
of singlets formed on alternating nearest-neighbor bonds 
(a pattern which can be realized in two different ways; 
hence the ground state is two- fold degenerate) Away 
from this special, exactly solvable point, there are fluc- 
tuations modifying the simple product state. Numerical 
exact diagonalization studies have shown that long-range 
dimerization survives down to gc ~ 0.241.—"— For g < gc 
the ground state exhibits critical spin and VBS correla- 
tions (like the standard Heisenberg chain with J2 = 0)r^ 
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FIG. 1: (Color online) A cylindrical, semi-periodic 2D square 
lattice with open edges (left and right sides) and periodic 
boundary conditions in the vertical direction (i.e., the open 
links at the top and bottom are connected to each other). 

At higher g, the simple dimer VBS order persists at least 
up to (7 sa 0.6, above which more complicated VBS or 
spiral spin states likely form.— 

The frustrated 2D square-lattice J1-J2 Heisenberg 
model (with nearest-neighbor couplings Ji and the J2- 
interactions connecting spins across the diagonals of each 
four-spin plaquette) also has a non-magnetic ground state 
in some window of coupling ratios 0.4 ^ g 0.6 (out- 
side of which the ground state is Neel antiferromagnetic 
for smaller g and exhibits stripe antiferromagnetic order 
for larger However, in this case it has been dif- 

ficult to determine the exact nature of the ground state. 
Many studies over the past two decades have suggested 
a VBS, with either columnar or plaquette (four-spin unit 
cell) order fi^"— but spin liquid ground states (which have 
no broken symmetries but may have topological order—) 
have also been proposedJ^*^ Very recently, results of 
DMRG calculations on cylindrical semi-periodic lattices 
(with open edges in one direction — see Fig. [T]) were used 
to argue more specifically that the ground state of the 
system for 0.41 < g < 0.62 is a Z2 spin liquidi^ A con- 
current calculation based on tensor-product states also 
claimed the absence of VBS order ^ 

A story similar to that of the J1-J2 Heisenberg model 
has played out in recent years in the case of the S = 
1/2 Heisenberg model with only nearest-neighbor inter- 
actions on the geometrically frustrated kagome lattice. 
Many calculations initially suggested a VBS ground state 
(in this case with a complex 12- or 36-site unit cell);^"— 
but the most recent DMRG studies support a Z2 spin 
liquid scenario^i^ (as also predicted in early analytical 
work2£). Here, as well, cylindrical lattices played a cru- 
cial role in obtaining the numerical data. 



B. Deconfined quantum-critical points and VBSs 
with emergent U(l) symmetry 

When the Neel order of a 2D antiferromagnet such as 
the 5* = 1/2 Heisenberg model is destroyed in a con- 
tinuous quantum phase transition, one scenario is that 
the putative spin-liquid state is immediately unstable to 
the formation of a VBS. This has been argued to lead 
to a "deconfined" quantum-critical point separating the 



Neel and VBS states i^ 1'^^ The phase transition is as- 
sociated with deconfinement of spinous. Being generi- 
cally continuous, due to subtle quantum interference ef- 
fectes, this type of transition violates the classical "Lan- 
dau rule" , according to which a transition between two 
ordered states breaking unrelated symmetries should be 
generically first-order. 

In the low-energy field-theory argued to describe the 
deconfined quantum-critical point (the 2+1 dimensional 
non-compact CP"'^ theory) (22 the VBS fiuctuations corre- 
spond to a U(l) gauge field to which spinous are coupled. 
There is a dangerously irrelevant operator (a quadrupled 
monopole operator) which reduces the U(l) symmetry to 
a four-fold {Z4) symmetry inside the ordered VBS state 
(in which the spinous become confined). On the square 
lattice, this corresponds to the four degenerate columnar 
VBS patterns. Close to the critical point, the Z4 sym- 
metry only becomes apparent beyond a length-scale A, 
which is larger than the standard correlation length ^ as- 
sociated with the magnitude of the order parameter. At 
distances below A there are angular fluctuations of the 
VBS order parameter (Dx,Dy), which in a system with 
X order, {\Dx\ > 0,Dy = 0), induces Dy order on length 
scales up to A, with this length diverging as A ~ ^1+° 
with a > 0. At distances much below A, the angle of the 
VBS order parameter fluctuates in an essentially U(l) 
isotropic manner. 

The deconfinement scenario appears to be realized in 
a class of "J-Q" modelsfii^"— in which the Heisen- 
berg exchange J is supplemented by certain multi-spin 
interactions — products of two or more two-spin singlet 
projectors acting on different spin pairs. These interac- 
tions lead to the formation of local correlated singlets, 
thereby reducing, and eventually destroying, the Neel 
order. Results of QMC calculations (which are not af- 
fected by sign problems in this case) are consistent with 
a single critical point separating the Neel state and a 
VBS. While some works suggested that the transition 
is weakly first-orderj^i^ the most recent studies point 
to a continuous transition with anomalously large scal- 
ing corrections 42r— Moreover, emergent U(l) symmetry 
has been explicitly observed in the VBS order parameter 
distributioujii^ By studying the U(l)-Z4 cross-over, the 
exponent a was estimated in Ref.|4l|to be a = 0.20±0.05. 



C. Stability of the spin liquid 

For the frustrated spin systems discussed above in 
Sec. II Al deconfined quantuni-criticality, i.e., a gapless 
spin liquid existing only at a singular point, is also an 
alternative to the transition out of the Neel state into 
an extended spin liquid phase. At the heart of this issue 
is the question of the stability of the spin liquid state4^ 
The deconfined quantum-criticality scenario implies that 
some spin liquids are generically unstable, at least un- 
der some commonly satisfied conditions, but stable spin 
liquids can also exist. 
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Recently Cano and Fendley succeeded in constructing 
a long-sought local (but complicated) Hamiltonian^ that 
is the parent Hamiltonian of the prototypical resonating 
valence-bond (RVB) spin liquid, i.e., the equal superposi- 
tion of all nearest-neighbor valence bond configurations 
(with the Marshall sign rule buih in)M-.^ This state, 
however, is a U(l) spin liquid with exponentially decaying 
spin correlations but critical VBS correlations j^Si^ and 
not the kind of fully gapped Z2 spin liquid proposed in the 
context of the frustrated models discussed above [but a 
U(l) spin liquid is also a possible ground state candidate 
of this model^] . As a consequence of its close relation- 
ship with the critical Rokhasr-Kivelson dimer model^^i^ 
one would expect this state to be generically unstable to 
perturbations of the Cano-Fendley Hamiltonian, leading 
to the formation of a VBS. Viewed from the perspective 
of a class of quantum states, the introduction of longer 
bonds either maintains the critical VBS?^ or leads to 
a Z2 spin liquidf^i^ but the hamiltonian for these ex- 
tended RVB states is not known. 

Stable Z2 spin liquids are known with Klein Hamilto- 
nians on particular decorated lattices;^ but the degree if 
stability of these states when moving away from the limit 
of high decoration is not known. The Kitaev honeycomb- 
lattice modelf^ which has a Z2 liquid state, can also be 
related to a model of SU(2) interacting spins on a deco- 
rated honeycomb lattice;^ However, there is still no rig- 
orously known example of a Z2 spin liquid ground state 
of a local SU(2) invariant Hamiltonian on one of the sim- 
ple standard 2D lattices (square, triangular, honeycomb, 
kagome, etc). This lack of a prototypical system under- 
lies the quest to find Z2 liquids in numerical studies of 
frustrated quantum spin Hamiltonians4i^— i^i^ Z2 spin 
liquid states have already been confirmed in QMC stud- 
ies of frustrated quantum XY models<^ 



D. Detection of spin liquids and VBS order 

It is highly non-trivial to unambiguously confirm 2D 
spin liquid states based on numerical calculations on rel- 
atively small lattices. The main difficulty here is to ex- 
clude weak VBS order (while the absence of magnetic 
order is easier to confirm, e.g., by demonstrating a non- 
zero spin gap). There is therefore much interest in find- 
ing positive signals for various spin liquid phases, e.g., 
using unique finite-size scaling properties of the entan- 
glement entropy ^^1^° Other signals related to the topo- 
logical aspects of spin liquids have also been proposed<^i^ 
However, regardless of what properties arc investigated, 
great care has to be taken in view of the small lattices 
accessible for systems with frustrated interactions. Due 
to sign problems, unbiased QMC studies of the ground 
states of these systems are essentially impossible^ (al- 
though some progress has been made here recently at el- 
evated temperatures^). Variational QMC methods can 
be used^i^ but are not reliable, because very different 
states can have almost the same energy. Exact diagonal- 



ization studies can reach 42 spinSfSs— while DMRG 
calculations now can reach hundreds of spins 1^ Tensor- 
product state methods (which can be regarded as gener- 
alizations of the matrix-product based^ DMRG scheme) 
can reach much larger sizes, but are complicated by the 
fact that extrapolations also have to be carried out in the 
bond dimension of the tcnsorsi^i^'^ In DMRG calcula- 
tions there is a similar issue with regards to the maximum 
number of states that can be kept, which is what limits 
the accessible system sizes (since that number of states 
in this case has to grow exponentially with the system 
size). 

E. Lattice shapes and boundaries 

As already mentioned, in DMRG studies it has become 
popular to use lattices in the form of cylinders with semi- 
periodic boundary conditions (with periodic boundaries 
along the long direction and open short edges), as illus- 
trated in Fig. [1] An aspect ratio L^/ Ly > 1 improves 
the convergence with the number of states kept, as com- 
pared to a fully periodic lattice with equal length in both 
directions (for a given total number of lattices sites) 1^ 
The better convergence with samples of this shape can 
be traced to the inherently ID nature of the DMRG pro- 
cedures and how the generated states can incorporate 
entanglement It has also been argued that cylindri- 
cal Lxl Ly > 1 samples, some times in combination with 
modifications of the boundaries (e.g., using field terms 
breaking some symmetry), have other favorable effects 
as well on the convergence of various order parameters 
as a function of the system sizcj^i^ 

In QMC studies of sign-problem free models, periodic 
L X L lattices are normally used. In cases where the 
couplings are spatially anisotropic, it has proved help- 
ful to use Lx X Ly lattices with ^ Ly^^^n^ while 
in other cases no particular advantages of such rectan- 
gular lattices were notedJ^ Open boundaries have been 
considered in QMC work primarily in cases where the 
perturbing effects of the edges are the actual targets of 
investigationiSiS In a previous QMC study of a VBS 
state it was also noted that open boundaries can be used 
to break the four-fold symmetry of the 2D VBS com- 
pletely and stabilize a unique VBS pattern, as an alter- 
native of studying VBS correlation functions in periodic 
lattices with no explicitly broken symmetriesi^ 

F. Purpose of the paper 

The main purpose of the present paper is to systemat- 
ically investigate the role of the lattice shape and bound- 
ary conditions on the finite-size scaling properties of the 
VBS order parameter. VBS states have in the past few 
years been conclusively demonstrated in several 2D J- 
Q modelsfi'^i^'2^1^ and also in ID chains (where the 
same kind of dimerization transition takes place as in 
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the frustrated J1-J2 chain) ^^1^^ and 3D systems Dif- 
ferent types of VBS patterns can be realized, depending 
on the arrangements of the singlet projectors on the lat- 
tice. These models have been studied with large-scale 
QMC simulations, mainly for the purpose of investigating 
the nature of the Neel-VBS transition.iii2,-ii>2^>26 Here 
the main focus will instead be on the VBS state itself 
(including its cross-over behavior close to criticality), us- 
ing the J-Q models to obtain generic bench-marks for 
finite-size scaling of this kind of order parameter. An ef- 
ficient approximation-free ground state projector QMC 
mcthod^i^ was used to obtain results for both strongly 
and weakly VBS ordered systems on square lattices with 
different shapes and boundaries. 

In order to make contact with the currently favored 
manner of applying the DMRG method,—"^ cylindrical 
systems with open edges in one direction will be stud- 
ied extensively. The convention adopted here is that the 
edges parallel to the y-axis are open, and periodic bound- 
ary conditions are applied in the other direction. Such 
an Lx X Ly lattice is illustrated in Fig. [T] In some cases 
the open edges will be modified to favor a certain VBS 
pattern, which is often also done in DMRG studiesj^i^ 
Fully periodic lattices will also be considered. Two as- 
pect ratios, L^/Ly ~ 1,2, will be considered for both 
the semipcriodic and fully periodic systems. The limit 
Lx ^ 00 will also be taken for small Ly. 

In addition to suggesting optimal approaches for ex- 
tracting the VBS order in the 2D thermodynamic limit, 
the results presented here will also show that the issue of 
excluding VBS order in a system with an unknown type 
of non-magnetic ground state may be more difficult than 
what has been anticipated so far. In particular, the ge- 
ometry of long cylinders can give misleading results. Not 
only can calculations on such systems completely miss 
2D VBS order (because the system is disordered with a 
short correlation length on the cylinders), but also the 
claimed positive signals of a 2D Z2 spin liquid^i^ can- 
not be trusted when used with cylinders of practically 
accessible widths (because they are essentially ID spin 
liquids although the state orders in the 2D limit). The 
emergent U(l) symmetry of the VBS state leads to in- 
teresting boundary effects, which are also studied here. 



G. Outline of the paper 

In Sec. |lT]the J-Q models are specified in detail, the 
correlation functions of interest are defined, the projec- 
tor QMC method is briefly outlined, and its convergence 
properties are discussed and illustrated with an exam- 
ple. Extrapolations of the inflnite-size value of the or- 
der parameter is discussed in Sec. IIIIl Results for the 
J-Q3 model at J = (the pure Q3 model), which has 
very robust columnar VBS order, is discussed first, in or- 
der to show how the different ways of extrapolating the 
order parameter to the thermodynamic limit agree well 
with each other. Results for three different lattice types 



are compared; periodic L x L and 2L x L systems as 
well as semi-periodic cylindrical 2L x L systems. The 
much weaker VBS ordering in the Q2 and J-Q2 mod- 
els is discussed next, using the same lattices as above. 
Here several subtle issues are pointed out that affect ex- 
trapolations to infinite size when the order is not strong, 
and, therefore, the length-scales ^ and A are large. The 
quantum-critical scaling form of the VBS order param- 
eter is also discussed, as a nearby critical point also in- 
fluences the finite-size behavior in systems off criticality. 
The vector aspects of the columnar VBS order parameter 
{Dj-, Dy) and the effects of its emergent U(l) symmetry 
are studied in detail in Sec. |lVl The evolution of the x- 
and y-components of the order parameter as a function 
of the distance from an open edge is studied, with and 
without symmetry-breaking modifications of the edge. In 
Sec.|V]the destruction of VBS order on cylinders is stud- 
ied in the limit 00 and Ly fixed. The most impor- 
tant results are summarized and their implications are 
discussed in Sec. |VT1 Here detailed comparisons with the 
recent DMRG results^ for the J1-J2 Heisenberg model 
are also made. Detection of the U(l)-Z4 symmetry of the 
VBS order parameter based on probability distributions 
P{Dx, Dy) generated in QMC calculations is discussed in 
Appendix |^ 

II. MODELS AND METHODS 

A. J-Q models 

A generic J-Q model is defined using products of sin- 
glet projectors C{i,j) on two sites, 

C(z,j) = i-S,.S,. (1) 

The standard Heisenberg model is just a sum of such 
singlet projectors over the interacting bonds (here 
nearest-neighbor sites on the square lattice), 

Hj^-jJ2Cii,j), (2) 

where the minus sign corresponds to antiferromagnetic 
interactions. A Q„ term consists of products of two or 
more (n) singlet projectors acting on different bonds; 

n 

a b=l 

Here a is a label corresponding to the lattice units within 
which the singlet projectors are arranged and b labels the 
bonds (spin pairs) on which the singlet projectors within 
these units act; i[a,b] and j[a,f'] above refer to the two 
sites connected by bond b in unit a. In the simplest kind 
of Q2 term on the square lattice, a denotes 2x2 pla- 
quettcs, with the two projectors within these plaquettes 
connecting spins either horizontally or vertically (i.e., for 
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FIG. 2: (Color online) Q2 and Qs terms on the square lattice. 
The bars of length one lattice constant indicate the locations 
of singlet projectors C{i,j) on site pairs The Hamiltonian 
contains all unique translations of these operators. 



a given 2x2 plaquette there are two labels a; one corre- 
sponding to horizontal and one to vertical bonds). This 
standard Q2 term will be considered here, along with a 
similar Qa term with the projectors arranged in columns. 
Both these cases are illustrated in Fig. [2l In general, the 
sum over projectors is such that the Hamiltonian does 
not break any of the symmetries of the lattice. 

The J-Qn model defined by the Hamiltonian H = 
Hj + Hq^ hosts a VBS ground state when Qn/ J is suf- 
ficiently large. In general, VBS formation is favored for 
a large enough number n of singlet projectors (with the 
minimum being typically n = 2 or n = 3 in two dimen- 
sions) if the arrangement of them is compatible with some 
symmetry-breaking pattern of strong and weak bond sin- 
glets. In this paper the pure Q2 and models without 
any J term will be studied primarily, but some results 
for J-Q2 systems with J/Q2 > will also be presented. 



when a Q term is present are straight-forward and have 
been discussed briefly in Ref. [lol ): First all the opera- 
tors in the J and Q terms are split into their diagonal 
and off-diagonal components in the basis of spin states 
jS'f , . . . , Sff) used. The diagonal operators can be moved 
around on the lattice as long as each operator is compati- 
ble with the spin state on which it acts (with only opera- 
tions on anti-parallel spins allowed) . The full set of oper- 
ators is sampled by changing the types of some operators 
from diagonal to off-diagonal, or vise versa, on the same 
lattice unit a, using an efficient loop algorithm.—"—). 

The ground state of a bipartite J-Q model (i.e., with 
each singlet projector connecting two spins on different 
sublattices) being guaranteed to be singlet, it is partic- 
ularly convenient to use a trial state expressed in the 
valence bond basis in the singlet sector. The conver- 
gence of {A)m to the true ground state expectation value 
(0|yl|0)m is then dictated by the gap to the second sin- 
glet. For a periodic lattice (or a semi-periodic cylinder), 
a transitional-invariant trial state also filters out excited 
states with non-zero momentum from the outset. Trans- 
lational invariance in the applicable lattice direction(s) is 
easily ensured by using an amplitude-product stated for 
1 5*0), i.e., a superposition written in terms of bipartite 
valence bond states 



(6) 



Here the sum includes all tilings of the TV-site lattice into 
N /2 bipartite two-spin singlets, i.e.. 



B. Projector QMC 

J-Q models with minus signs as in Eqs. ^ and ([3]) do 
not have QMC sign problems and can be studied with 
very efficient QMC loop algorithms. Here the ground 
state projector method developed in Ref. Ull is used. It 
is based on applying a high power of the Hamiltonian to 
a "trial" state l^'o). 



(4) 



where (— iJ)™ is written as a sum over all possible strings 
of the individual J and Q terms in ^ and ([3]). Denot- 
ing such a string of singlet projectors by P„i{i), with i 
formally indexing the different strings, an operator ex- 
pectation value is written as 



(5) 



where P^ij) is the string Pm{j) in reverse order. 

The QMC method implements importance sampling 
of the operator strings P^{j)Pm{i), which is done in two 
steps, as described in detail in Ref. [8l| in the case of the 
Hcisenbcrg model (and the modifications of the scheme 



N/2 



(7) 



where (i, j) = (| tiij) — I iitj))/V^ with i and j sites on 
sublattice A and B, respectively, and the weight Cy of a 
given tiling v into singlets depends only on the "shapes" 
1 = {Ix, ly) of the bonds in \v); 



(8) 



where ni is the number of bonds of shape 1. 

Amplitude-product states arc very easy to sample in 
the course of the projection according to (jlj, as also de- 
scribed in Ref. |8l|. The detailed form of the amplitude 
h(l) is not crucial when the state is used as a trial state. 
Variationally optimized amplitudes lead to faster conver- 
gence with the power m, but even without optimizing the 
convergence properties are good<^ In the work reported 
here, amplitudes decaying with the bond length I as 
were used (in which case the trial state itself has Neel or- 
der, but this is very quickly destroyed by the projection 
procedure in a VBS state). 
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C. Correlation functions 

In order to characterize the ground state, the spin (s) 
and dimer (d) correlation functions are computed. These 
are defined in the standard way as 

Cs{v,j) = (S(rO-S(r,)), (9) 
C^da(r,j) = (B„(r,)S„(r,)), (10) 

where r^j = r,; — Vj is the spatial separation of the opera- 
tors and Ba, a ~ x,y, is the dimer operator on nearest- 
neighbor bonds oriented in the a direction, e.g., for a ~ x 

B£(r) -S(r).S(r + i). (11) 

One can also define cross-correlations {Bx{vi)By(Tj)) but 
they will not be needed here. The correlation functions 
can be easily computed using loop estimators based on 
the transition graphs generated when the sampled va- 
lence bond states in the ket and bra states of Eq. ([5]) 
are propagated by the string of singlet projectors. The 
estimators are discussed in detail in Refs. [73 andlssl. 

Columnar and plaquette VBS states can both be de- 
tected by the columnar VBS order parameter, which 
when averaged over the whole lattice of = L^Ly sites 
can be defined by the operators 

= (12) 

= ^E^s(-^>y)(-iF- (13) 

In a columnar state with the lattice rotational symmetry 
completely broken, either or Dy has a non-zero ex- 
pectation value, while in a plaquette state they are both 
non-zero and equal. The J-Q models studied here host 
only columnar VBS states. However, as we will be dis- 
cussed below, in a columnar state on a finite lattice one 
can have both non-zero (Dx) and (Dy), due to boundary 
and shape effects. 

In periodic and semi-periodic systems where the de- 
generacy of the possible VBS patterns is not broken, one 
can only detect the VBS with the corresponding correla- 
tion functions, e.g., the squares of the order parameters 
defined above. In particular, it is useful to consider the 
total squared order parameter, 

D^ = Dl + Dl. (14) 

The magnitude of the order parameter in a correspond- 
ing symmetry-broken state is _D = (_D^)^/^ (which can 
be taken as a definition of the value D of the order pa- 
rameter) . In non-square samples it is also illuminating to 
investigate the components (D^) and {Dy) individually, 
to see how the lattice shape (and boundaries) affect the 
symmetry breaking. As will be demonstrated in the fol- 
lowing sections, this issue is, in fact, of key importance 
for interpreting numerical results for non-square samples. 



In the cylindrical semi-periodic systems it is useful to 
define the order parameter in such a way that the per- 
turbing effects of the open edges are partially eliminated. 
As in Ref. 0, for such systems with > Ly the summa- 
tions in and p3| will normally be taken over only 
the central sites within a square of size Ly x Ly. 

In cases when the lattice coordinates (x, y) are needed 
explicitly in the further discussion of correlation func- 
tions in the later sections, the numbering convention will 
be a; G {0, . . . , — 1} and y G {0, . . . ,Ly — 1}. 

D. Convergence tests 

To examine the convergence properties of the projector 
method, the state Q after m operations with H can be 
written in terms of eigenstates \n) oi H as 

IV-™) (15) 

n 

where c„ are the expansion coefficients of the trial state in 
the energy basis. Assuming that the ground state energy 
Eq is the largest in magnitude, |i?o| > \En\, V n > 0, 
which is the case for sure with a Hamiltonian expressed 
using the singlet projectors ^ and the signs as in ^ and 
an expectation value of an operator A not commuting 
with the Hamiltonian can be expanded as 

(A),„ = (0|A|0)+2(l|A|0)^(^-^j +.... (16) 

Here |1) is the first excited state in the symmetry sector 
considered, which with an amplitude product state obey- 
ing all applicable lattice symmetries is a singlet that is 
fully symmetric with respect to all the symmetry oper- 
ations (translations, reflections, and rotations of the lat- 
tice) . With the gap A = Ei — Eq and a large projection- 
power m, Eq. (jl6|) can be written as 

(A),„ = (0|A|0)+cxexp(^-^A^, (17) 

where cq is the ground state energy per site, eo — Eq/N, 
and c is a constant. In order to achieve good conver- 
gence, one should therefore use a size-normalized projec- 
tion power m/N ^ 1/A. 

The gapped VBS state being of primary interest here, 
A/eo approaches a non-zero constant as the system size 
increases. One may then expect good convergence prop- 
erties with an essentially size independent m/N . How- 
ever, for system sizes accessible in practice, the gap still 
typically decreases significantly with the system size. In 
addition, the density of states above the gap increases 
as well. As a consequence, m/N has to be increased 
with the system size to ensure good convergence. Since 
the number of operations required for one full sweep of 
Monte Carlo updates of a configuration in the projector 
method is of order the computation time in practice 
grows faster than N . 
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FIG. 3: (Color online) The VBS order parameter as a function 
of the projection power m (normalized by the system size A'^) 
in simulations of the Q3 model on a periodic 32 x 32 lattice. 
The inset shows an exponential fit of the form (|17|l . 

All results presented here were tested for convergence 
by carrying out several calculations with different pro- 
jection powers m/N oc L (with L = m.ax[Lx, Ly\ for 
non-square lattices) and making sure that there is no 
remaining detectable dependence on m. An example of 
a detailed convergence test is shown in Fig. [31 Typically, 
m/N = L/2 was sufficient to ensure good convergence. 
In principle the singlet-singlet gap can be extracted by 
fitting the exponential form ([T7)) to data such as those in 
Fig. [3] (as shown in the inset), but such gaps will not be 
studied here. 



III. EXTRAPOLATION OF VBS ORDER 

Previous ground-state and finite-temperature QMC 
studies have confirmed that both the J-Q2 and J-Qs 
models, with the singlet projectors arranged as in Fig. [2j 
have VBS-ordered ground state for large Qn/JJ^'^'^ 
The maximal order parameter obtains for J ~ Q (pure 
Qn models) and, naturally, the order is more robust in 
the Qs model. The previous studies were mainly con- 
cerned with the critical and near-critical aspects of the 
Neel and VBS order parameters — the critical exponents 
as well as the emergent U(l) symmetry seen in the VBS 
order parameter (Dx,Dy). 

In this section some important aspects of the VBS or- 
der parameter will be discussed first, in particular the 
expected consequences of its emergent U(l) symmetry. 
Then, turning to numerical results, the magnitude of 
the VBS order parameter of the pure Qs model will be 
extracted first, to illustrate the convergence as a func- 
tion of the lattice size for several cases of lattice shapes 
and boundary conditions. The J-Q2 model, including 
the pure Q2 case, is then considered in order to inves- 
tigate potential problems arising when the VBS order 
is weaker. The quantum-critical scaling will also be dis- 
cussed briefly, as it is directly related to the extrapolation 



problems when the VBS can be considered near-critical. 



A. Nature of the VBS order parameter 

Note first that the maximal columnar VBS order pa- 
rameter is obtained for the state with no fiuctuations in 
the valence bond basis — the state with nearest-neighbor 
singlets on all bonds of every second column. If the sin- 
glets arc oriented in the x direction, then the order pa- 
rameter components defined in (jl2[) and (|13p have the 
expectation values (Dx) =3/8 (up to an arbitrary sign) 
and (Dy) = 0. If the symmetry is not broken and the 
ground state is an equal superposition of the four de- 
generate valence-bond states with horizontal and verti- 
cal bonds, then the expectation value of the squared VBS 
order parameter HI]) is {D'^) = + {Dy^ = 9/64 in 

the limit of an infinitely large system. For finite systems 
there are corrections to this value, however, which are 
related to the non-orthogonality and over-completeness 
of valence bond statesi ^°'^^ 

The emergent U(l) symmetry property of the VBS or- 
der parameteri^"— and its related length-scale A will 
be of importance in order to understand many of the re- 
sults to be discussed here and in the later sections. For 
i ^ A, the order parameter {Dx,Dy) on an L x L lat- 
tice behaves essentially as an isotropic 2D vector, while 
for X ^ A the order parameter locks to one of the four 
angles mr/2. This is further discussed in Appendix [XI 
Here, for 7^ Ly lattices, with or without open edges, 
the U(l)-Z4 cross-over will manifest itself also in how 
(on what length scale) the 90° rotational symmetry of 
the VBS order parameter is broken on a lattice which 
does not have this symmetry. 

It should be noted that symmetry cross-overs such 
as the U(l)-Z4 case discussed here also occur in many 
classical systems with dangerously irrelevant perturba- 
tions (i.e., ones that do not change the universality 
class of a phase transitions but reduce the degeneracy 
of the ordered state), e.g., the 3D XY- model with a q- 
fold symmetry-breaking field of the form cos{q9i) (with 
q > 4)4^1^ There are several numerical studies of the 
scaling dimension of such a dangerously irrelevant per- 
turbation and the nature of the cross-over and its length- 
scale A^-^ 



B. Strong VBS order in the pure Q3 model 

Fig. [H shows the size dependence of {D'^) of the Q3 
model computed on periodic L x L and 2L x L lattices. 
For the latter systems the individual expectation values 
{Dx)'^ and (Dy)'^ are also shown (while these are of course 
both equal to (D^)/2 for the L x L lattices). Here the 
convergence of {D'^) to a non-zero value when L ^ 00 
is apparent for both types of lattices. It is interesting to 
note that both the x and y components are nonzero on 
the 2L X L lattices for small L, but for larger systems 
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FIG. 4; (Color online) Size dependence of the squared order 
parameter and its x and y components of the Qa model com- 
puted on periodic Lx L and 2Lx L lattices. The curves pass- 
ing through the (D^) data are second-order polynomial fits 
(excluding the systems for which this form cannot be used). 
The error bars are much smaller than the plotting symbols 
(typically the standard deviation is ~ 2 x 10^^). 

the symmetry is completely broken, eventually leading 
to (Dl) 0, (Dy) {D"^)- Thus, on the non-square 
periodic lattices the columnar state with the bonds ori- 
ented parallel to the shorter lattice direction (here Ly) is 
energetically favored. This remains true also for larger 
aspect ratios Lx/ Ly. 

The cross-over from partially broken to fully broken x- 
y rotation symmetry, which in Fig. 2] takes place for the 
2LxL systems for L 20, should be related to the emer- 
gent U{1) symmetry of the VBS order parameter . ^'^^i'^^ 
As discussed in Appendix \^ for the Q3 model no per- 
fect U(l) symmetry can be detected on periodic L x L 
lattices (since the length-scale A is very short), but for a 
wide range of sizes the system is in a cross-over regime be- 
tween U(l) and Z4 symmetry. The range of L over which 
the cross-over to a purely y-ordered VBS takes place in 
Fig. |4]is roughly where all traces of U(l) symmetry van- 
ish on the L X L lattices (as discussed in Appendix [X^ . 

Turning now to the quantitative behavior of the to- 
tal order parameter for the largest systems in Fig. |31 as 
expected the order parameters for both lattice types ex- 
trapolate to the same value in the thermodynamic limit. 
Fits of the data for the largest systems to second-order 
polynomials are shown. Note, however, that this form is 
strictly not correct. For a discrete broken symmetry one 
would expect the asymptotic finite-size corrections to be 
exponentially decreasing with increasing system size. It 
is not easy to reach sufficiently large systems to observe 
this behavior, however. The second-order fits look rea- 
sonably good on the scale of the plot, but in fact they are 
not of high quality statistically when 6-8 data points are 
included. Including higher powers helps somewhat, but 
this can lead to fitted forms that do not behave monoton- 
ically as 1/L 0. Such problems with the polynomial 




1/L 

FIG. 5; (Color online) Size dependence of the squared order 
parameters of the Qs model on cylindrical 2Lx L lattices (us- 
ing the central L x L square for computing the expectation 
values) . The smooth curves are second-order polynomials fit- 
ted to the (D^) data for several of the largest system sizes. 



fits reflect a cross-over to the eventual exponentially rapid 
convergence. Using second-order fits for the largest few 
system sizes still should result in a reasonably accurate 
extrapolated order parameter. Normally such an extrap- 
olation should give a lower bound on the actual value, but 
this cannot be guaranteed in the presence of statistical 
errors. In the case considered here, the results for L x L 
and 2LxL extrapolate to 0.0691 and 0.0684, respectively, 
with the fits shown in Fig. Because of the issues with 
the, strictly speaking, wrong form of the fitting func- 
tion, it is not meaningful to compute error bars on these 
numbers — the purely statistical errors are smaller than 
the variations among fits with different polynomials and 
number of data points included. For the purposes of the 
investigations in this paper, the issue of statistical errors 
is only of minor importance, however (while the system- 
atical errors due to a wrong fitting form are important). 

Data for cylindrical 2L x L systems are shown Fig. [S] 
Here the order parameters arc computed on the central 
L X L square. In sharp contrast to the fully periodic 
2LxL systems, here it is the x component of the order pa- 
rameter that survives in the thermodynamic limit. Thus, 
the open edges along the y direction favor the bonds or- 
dering perpendicularly to them, and this effect wins over 
the competing effect, demonstrated in Fig. 01 of the as- 
pect ratio Lx/Ly > 1 favoring bonds ordering in the y 
direction. Quadratic fits to (D^) and (D^) for a few of 
the largest system sizes extrapolate to 0.0685 and 0.0694, 
respectively, in good agreement with the results for the 
periodic systems. 

As a consequence of the open boundaries inducing an 
x-oriented VBS, the ordering pattern in this case is non- 
dcgencratc. Therefore, the unsquared order parameter 
{Dx) is non-zero and should, in the thermodynamic limit, 
take a value agreeing with the squared order parameters 
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FIG. 6: (Color online) Local columnar x order parameter (|18|) 
of the Q-i model computed at the center of a 21/ x L cylinder. 
The smooth curve is of the exponential form (|17|l and extrap- 
olates to 0.264. The inset shows the location dependent bond 
correlation function {Bx{x)) for a 32 x 16 system. 

extracted above; {D^) (Z?^)^/^. The expectation value 
of the nearest-neighbor spin correlator ([TT|) indeed oscil- 
lates considerably as a function of the location along the 
X direction, as shown in the inset of Fig.lHlfor the 32 x 16 
cylinder. The dimer order is clearly the strongest at the 
edges but remains large also in the interior of the system. 

A local VBS order parameter for a system with bonds 
ordered along the x axis can be defined as 

D^ix) = {B,{x,y))-\{B,{x-l,y) + B^{x + l,y)), (18) 

which is independent of y on the semi-periodic cylindrical 
lattices (and can be averaged over y in the QMC calcu- 
lations). This quantity at the central column is shown 
as a function of the inverse system size in the main plot 
of Fig. [51 Here an asymptotic exponentially fast conver- 
gence can be seen clearly, which is illustrated with a fit to 
the form ([T7]). This fit is of good statistical quality and 
extrapolates to 0.264, in good agreement with the values 
for {jfi)^!"^ obtained above. The magnitude of the order 
parameter of the model is, thus, 70% of the largest 
possible value (3/8) for a columnar VBS. 

C. Reduced order in the Q2 and J-Q2 models 

In the pure Q2 model, the VBS order is considerably 
weaker than in the Qj, model. The first study of this 
model gave the order parameter D w 0.070, or about 
20% of the maximal value, based on extrapolations of 
L X L results for L < 32^ While this order may still be 
regarded as quite strong, problems with extrapolating it 
correctly based on small to moderate lattice sizes already 
start to become apparent. 

Fig. [7] shows results for periodic L x L systems with 
4 < L < 72. A 5th-ordcr polynomial can be fitted very 
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FIG. 7: (Color online) Size dependence of the squared VBS 
order parameter of the pure Q2 model on periodic L x L 
lattices. The solid black curve in the main graph shows a 
fit of the L < 12 data to a second-order polynomial (which 
extrapolates to an unphysical negative value when L — > 00). 
The solid red curve shows a 5th-order polynomial fit to all 
the data, while the dashed black curve shows a quadratic fit 
to only the L > 20 data. The inset shows the behavior for 
the largest systems on a more detailed scale. 

well to all these data and extrapolates to 0.0063, about 
10% lower than the previous result. However, if only 
i > 20 data are used, a second-order polynomial is suf- 
ficient and the extrapolated value is significantly lower; 
(D^) = 0.0058. This illustrates the fact that polynomial 
fits based on small systems are not very reliable, because 
of the eventual exponential convergence (which is not yet 
fully apparent for the system sizes accessible). The re- 
sulting relative uncertainties are much larger than in the 
strongly ordered Qs model. The extrapolated value de- 
pends significantly on what system sizes are included in 
the fit and the order of the polynomial used. For the 
system sizes studied here, a pure exponential form does 
not yet work. 

An important aspect of the finite-size scaling behavior 
in the Q2 model is that the data for small to moderate 
lattices do not even clearly point to an ordered ground 
state. Fig. [7] also shows a second-order fit to only the 
i < 12 data points. The fit is statistically sound, but 
extrapolates to a negative value. Without access to larger 
system sizes it is not possible, using fitting procedures 
like this in 1/L, to determine whether the ground state 
of the infinite 2D lattice is ordered or disordered. At least 
i = 20 is needed with 1/L extrapolations to definitely 
conclude that the ground state is ordered. It can be noted 
that an asymptotic oc 1/L^ behavior is expected if there 
is no long-range order, but this form should apply only 
for L much larger than the correlation length. Note that 
the correlation length itself is also not easy to extract 
from the correlation functions unless L ^ ^ (which is 
not the case here). 

Fig.[5]shows results for periodic 2Lx L lattices. Using 



10 




1/L 

FIG. 8: (Color online) Size dependence of the squared total 
VBS order parameter {D'^) and its individual x and y compo- 
nents of the pure Q2 model on fully periodic 2L x L lattices. 
The curve displayed for 1/L < 0.125 is a 4th order polynomial 
fit to the (D^) data for L > 16. 



polynomials to reliably extrapolate results to the infinite 
size limit is again difficult. An example using a 4th order 
polynomial with data for L > 16 is shown which extrap- 
olates to {D^) = 0.0063. Here it is again clear that the 
polynomial is not the correct form, because the fitted 
curve deviates significantly for the smaller systems not 
included in the fit. 

The behavior of the individual x and y components in 
Fig. IS] appears to be qualitatively different from that ob- 
served in the Q3 model (Fig. 2]). In the more strongly or- 
dered Q3 model the y component is always significantly 
larger than the x component, and for large systems it 
completely dominates (the x component vanishing). In 
the Q2 model, the length-scale A of the cross-over from 
U(l) to Z4 symmetry is much larger, and the dimer or- 
der parameter acts as an essentially isotropic vector even 
for the largest lattices considered here. A cross-over to 
a behavior where the x component vanishes (as in the 
Qs model) should take place for larger system sizes, but, 
according to the analysis for L x L lattices in Appendix 
\X[ the cross-over length is beyond what can currently be 
studied with QMC calculations, with there being only 
weak signals of a columnar state. Since the two compo- 
nents arc almost equal in magnitude in Fig. |51 not know- 
ing about the peculiar finite-size effects due to emergent 
U(l) symmetry one may draw the erroneous conclusion 
from these data of the system being a plaquette VBS. 

It is interesting to note in Fig. [8] that the emergent 
x-y symmetry is not manifested yet for the smallest sys- 
tems. This reflects that fact that the continuous angular 
nature of the VBS order parameter only appears upon 
coarse-graining and L < 10 is not sufficiently large for 
representing a continuous VBS angle. The two cross- 
over length-scales, into and out of an U(l) symmetric 
order parameter, have been investigated in detail in clas- 
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FIG. 9: (Color online) Size dependence of the squared total 
order parameter (D^) and its x and y components, computed 
for the Q2 model on cylindrical 2Lx L lattices (including only 
the spins on the central L x L square in the definition of the 
order parameters). The curves are polynomial fits. In the 
case of Dy, no constant term was included. The inset shows 
the data for large systems on a more detailed scale. 

sical systems (clock models) exhibiting emergent U(l) 
symmetry)^ 

As in the Q3 model, on the open-edge cylinders with 
Lx > Ly the favored VBS ordering pattern is that with 
the bonds primarily in the x-direction. Fig. [5] shows re- 
sults for 2Lx L cylinders. Here the effect of the edges to 
strongly favor x ordering overcomes the tendency to U(l) 
symmetry, and there is never any size range for which the 
X and y components are almost equal. Also here the be- 
havior of both components for small lattices exhibit a 
naive extrapolation to a negative order parameter. For 
larger lattices (D^) crosses over to a form extrapolating 
clearly to a non-zero value, while the y component ex- 
trapolates to zero. A 4th-order polynomial fit to all the 
x-component data gives (-D^.) = 0.0047. This is signifi- 
cantly lower than the value quoted above for the exam- 
ples of extrapolations oi L x L data. However, the ex- 
trapolation is again sensitive to the lattice sizes included 
and the form of the fitting function used. 

It is also useful to examine the long-distance VBS cor- 
relation function, which should contain less finite-size 
corrections to the infinite-size order parameter than the 
sums over all correlations. The squared order param- 
eter p4|) contains significant non-asymptotic contribu- 
tions from short distances. Using the real-space dimer 
correlation function defined in Eq. (fT0|) . the staggered 
part in the case of the x component (and an analogous 
form for the y component) can be extracted according to 

QJ.T,y) = iQ,(x,y) (19) 

-jiCdAx - 1, y) + Cdx{x + 1, y)], 

where a factor 1/2 has been included in order for 
^ oo,y) {Dl), with At defined in Eq. 
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FIG. 10: (Color online) Staggered component, Eq. (|19|l of the 
long-distance dimer correlations in the Q2 model on periodic 
LxL and 2LxL lattices. Here r is the longest distance on the 
lattices; r — {Lx/2, Ly/2). The curve through the LxL data 
is a high-order polynomial fit. The inset shows the 2L x L 
data on a more detailed scale. 

in the thermodynamic limit. Fig. [10] shows results for 
the longest distance on periodic LxL and 2L x L lat- 
tices. For the LxL systems the sum of the x and y 
components is shown, along with a high-order polynomial 
fit that extrapolates to the infinite size order parameter 
= 0.0061. This extrapolation should be reasonably 
reliable, because the data for the largest systems flatten 
out clearly, reflecting the asymptotic exponential con- 
vergence (unlike the integrated quantity (D^) in Fig. [3 
where no flattening-out is yet seen). For the 2L x L sys- 
tem no reliable extrapolation is possible, because both 
components exhibit non-monotonic behavior. The sum 
of the X and y correlations for large L is nevertheless 
very close to the LxL results. 

It can also be noted in Fig. [10] that the individual com- 
ponents of the correlation function at long distance show 
somewhat less prominent x-y symmetry than the inte- 
grated correlators in Fig. H] although they are both still 
roughly equally large. Again, in the thermodynamic limit 
one of the components, likely the x component, will have 
to turn down and vanish, as the L^ > Ly geometry favors 
ordering in the y direction. 

In the previous section, it appeared that the most re- 
liable way to extract the order parameter in the thermo- 
dynamic limit is to exploit the symmetry-breaking open 
edges, using (D^(x)) defined in Eq. ([TS]) . Fig. ITT] shows 
such results for open-edge cylinders of size L x L as well 
as 2L X L. Here the induced x order appears to ex- 
trapolate to a value below the one obtained in Fig. [10] 
based on the long-distance correlation function — for the 
2L X L systems {Dx{L — 1)) is almost size indepen- 
dent for the largest systems, and one might hence con- 
clude that it has converged. The square of this value is 
(Dx)'^ sa 0.073^ « 0.053, which seems too low compared 
to the results in Fig. [TUl 




FIG. 11: (Color online) Main graph: Open-edge induced or- 
der parameter of the Q2 model at the center of cylindrical 
LxL and 2L x L lattices. The horizontal line corresponds 
to the value of the infinite-size order parameter from the ex- 
trapolation in Fig. 1101 



The reason for this apparent inconsistency should 
again be related to the emergent U(l) symmetry of the 
VBS order parameter: In addition to the x compo- 
nent of the order parameter induced by the open edges, 
there still remains, for the accessible lattice sizes, a non- 
negligible y component. This component is not locked- 
in by symmetry-breaking boundaries, however, but aver- 
ages to zero if measured without first taking the square of 
its operator. The existence of a non-negligible fluctuating 
y component nevertheless reduces the induced (Dx) from 
the full value, which should satisfy = (D^)^ + {Dl) 
for large systems. It is only when the system size exceeds 
the U{\) length scale A that one can expect the full or- 
der parameter to condense into the component (Dx), and 
this length scale cannot at present be reached for the Q2 
model. This shows again that the problem of extracting 
the VBS order parameter in the thermodynamic limit is 
a very delicate one. 

The examples shown here demonstrate that, when the 
VBS order is relatively weak (the length scale A is large), 
it is important to look at the full order parameter, in- 
cluding both the X and y components. The long-distance 
correlation function (Fig. ITU)) on L x L periodic lattices 
seems to be the fastest converging quantity, and it is in 
most cases best to use LxL lattices for extrapolations. 

When turning on the Heisenberg exchange J, the 
VBS order of the J-Q2 model is reduced and vanishes 
when J/Q2 ~ 0.0454^ Here two cases are considered, 
J/Q2 = 0.03 and 0.10, with the latter corresponding to 
a near-critical Neel state. Fig. [T^] shows results for the 
total squared VBS order parameter and the staggered 
part of the dimer correlation function (|19p averaged over 
the X and y directions. With (I?^) graphed versus 1/L 
it is again difficult to extrapolate to infinite size based 
on small lattices. Here the lattices are nevertheless suf- 
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FIG. 12: (Color online) The staggered part, Eq. (O, of the 
long-distance correlation function (at rmax ~ V^L) and the 
total dimer order parameter for the J-Q2 model at J/Q2 ~ 
0.03 and 0.10 on periodic L x L lattices. 



ficicntly large for it to be apparent that the system at 
J/Q2 = 0.03 is VBS ordered, while for J/Q2 = 0.1 the 
decay is much more rapid and consistent with no VBS or- 
der. The corresponding long-distance correlations show 
these behaviors much more clearly, with the J/Q2 = 0.03 
data exhibiting the expected exponentially fast conver- 
gence to a non-zero value for the largest sizes. Still, if 
data only for L up to w 10 were available, it would not be 
possible to unambiguously confirm the presence of long- 
range VBS order, even though the order parameter here 
is still above 10% of the maximum value. 

Note that the long-distance correlation function de- 
cays exponentially as a function of 1/L in a non-VBS 
state, i.e., much faster than the 1/L^ behavior of the total 
squared order parameter. It is therefore also much easier 
to confirm the absence of long-range order by studying 
the long-distance correlations. 



D. Quantum-critical scaling 

Ultimately, the difficulties in extrapolating the VBS or- 
der parameter to infinite size based on small systems will 
in many cases be related to critical scaling in the proxim- 
ity of a quantum-critical point (or "pseudo-critical" scal- 
ing in cases where the transition out of the VBS state 
is weakly first-order). A small system exhibits quan- 
tum criticality also slightly away from the critical point. 
Hence, data for a series of lattices may appear to extrapo- 
late to a disordered state, even though the infinitely large 
2D system is on the VBS side of a quantum phase transi- 
tion. According to conventional finite-size scaling theory, 
the window around the critical point within which a sys- 
tem of linear size L exhibits scaling is proportional to 
L'^/" , where v is the exponent governing the divergence 
of the correlation length. Depending on the prcfactor, 
this window may be sizable for practically reachable lat- 
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FIG. 13: (Color online) Size dependence of the VBS (top) and 
Neel (bottom) order parameters of the J-Q2 model at four 
different coupling ratios. The point J/Q2 = 0.0447 should 
be very close to the quantum-critical value according to the 
scaling analysis of the spin stiffness carried out in Ref. |4^ . 
The straight lines fitted through the J/Q2 = 0.0447 data (for 
system sizes L > 32) have slope —1.27 in both cases. 



tice sizes. As will be shown next, this is one reason why 
fits to small-lattice data can give misleading results, e.g., 
in the case of Q2-inodel results in Fig. [T] 

In addition to illustrating the near-critical VBS, the 
scaling of the Neel order parameter will also be briefly 
discussed here. According to past studies, both the J- 
Q2 and J-Qs models are strong candidate o^'^^ for the 
deconfined quantum-criticality scenario^ according to 
which both order parameters should be critical exactly 
at the same point. Results for the J-Q2 model will be 
discussed here. 

While all numerical results so far are consistent with 
a single Neel- VBS transition point, it has proved re- 
markably difficult to determine the location {J/Q2)c of 
this transition precisely. The most recent QMC stud- 
ies point to a continuous transition with unusually large 
scaling corrections in the quantities normally used to ex- 
tract the critical point, e.g., the spin stiffness and Binder 
cumulants42r— These corrections have made it difficult 
to reliably extrapolate the critical coupling ratio {J/Q2)c 
to infinite size. By using a logarithmic scaling correc- 
tion to the spin stiffness (which was not predicted in the 
original field-theory description of deconfined quantum- 
critical points but may appear with a modified actioi*^), 
{J/Q2)c = 0.0447 ± 0.0002 was obtained in Ref . [H. Us- 
ing a conventional correction oc L~", with small oj and a 
large prcfactor (which potentially could be a consequence 
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of the dangerously irrelevant operator responsible for the 
Z4 symmetric VBS), gives a similar result. 

In Fig. [13] the two order parameters are graphed versus 
the system size on log-log scales for coupling ratios close 
to the critical value. The Neel order parameter (Af ^) (the 
squared sublattice magnetization) is the size-normalized 
(tt, tt) Fourier transform of the spin correlation function 
([5]) . Both order parameters indeed exhibit critical scaling 
at {J/Q2) = 0.0447. For other couplings the curves fan 
out in the way typical for critical points. 

Interestingly, at J/Q2 = 0.0447 both order parameters 
scale as L~(i+'') with 77 « 0.27 (with a purely statistical 
error bar of about 0.01) when L < 32 systems are used 
in the fits. For the sublattice magnetization, this expo- 
nent is slightly smaller than in previous works , ^'^1^^ while 
the VBS exponent is somewhat larger than in Refs. [l| 
and nil If these exponents are truly exactly the same, it 
would imply a duality of the effective low-energy field the- 
ory that had not been anticipated,— but further detailed 
work, using larger system sizes and studying several cou- 
pling ratios in the neighborhood of J/Q2 = 0.0447, will 
be required before such a claim can be made (and it could 
also be a coincidence that the two exponents are almost 
equal). Note also that the value of quoted here may 
also still be affected by sub-leading scaling corrections. 

For coupling ratios larger than the critical value, in 
Fig. [Ql exemplified by J/Q2 = 0.1, the VBS order pa- 
rameter turns downward, reflecting the faster decay to 
zero. Asymptotically, in the Neel state the decay should 
follow the 1 /L"^ form, but this can only be observed when 
the lattice size exceeds the correlation length (which is 
very large this close to the critical point). The sublattice 
magnetization turns upward, reflecting an extrapolation 
to a non-zero value. For smaller J IQ-i-, here and 0.03, 
the behavior is the opposite, reflecting a VBS state with 
no coexisting VBS order. 

For the present purpose of detecting VBS order, an 
important aspect of the critical scaling is that, once a 
critical point has been identified, upward deviations from 
the power-law scaling, as seen in Fig.[T3]at J/Q2 ~Q and 
0.03, still can demonstrate an ordered state when moving 
away from criticality. It may be easier, in many models, 
to establish a critical point (or a first-order transition) 
than to accurately extrapolate the infinite-size value of 
the order parameter in a state with significant fluctua- 
tions (an order parameter significantly smaller than its 
maximum possible value). Based on the knowledge of 
the existence of a phase transition, it may be possible to 
establish long-range order even in the presence of strong 
quantum fluctuations. This will be the case especially in 
calculations limited to much smaller systems. 



IV. BOUNDARY SYMMETRY BREAKING 

One interesting aspect of the results presented in the 
previous section, exemplified in Figs. H] and El is that the 
boundary conditions dictate which of the order parame- 



FIG. 14: (Color online) A cylindrical lattice with modified 
open edges favoring columnar VBS order with vertical bonds. 
The thick vertical bars represent Q2 terms excluded from the 
summation in the Hamiltonian, Eq. Q. 

ter components, {D"^) or {D"^), is the one surviving in the 
thermodynamic limit. For 90° rotationally-symmetric 
periodic LxL lattices both order parameters are of course 
equal by symmetry (and spontaneous symmetry break- 
ing in the thermodynamic limit will randomly select one 
of the directions), but in other cases only one of them 
should survive in the thermodynamic limit (i.e., the lat- 
tice shape acts like a symmetry-breaking field). Exactly 
how the symmetry is broken should be model dependent, 
and also dependent on fine details of the boundary condi- 
tions. Note that there are no "neutral" boundaries for a 
VBS, i.e., any boundary conditions should favor one com- 
ponent of the order parameter above the other (expect 
perhaps for some unusual fine-tuned boundaries with ad- 
justable couplings). 

Here the Q2 and Q3 models will be used to illustrate 
the complexity of the boundary issues further, with direct 
measurements of the order-parameter components {D^) 
and {Dy) in systems where the edges break either the a;- 
translational symmetry or both the x and y symmetries. 
The boundary effects are particularly interesting in view 
of the emergent U(l) symmetry, due to which both order 
parameter components can survive up to large system 
sizes, as already shown in Scc. lIIIBj in the case of periodic 
systems. Here the ability of boundaries to twist the local 
order parameter in the {Dx^ Dy) plane will be studied. 

Two types of 2L x L cylindrical lattices will be used. In 
addition to the case discussed so far, where the y-oriented 
edges are open and uniform, a modified boundary that 
breaks the translational symmetry in the y direction will 
also be studied. The modification acts as a field induc- 
ing Dy order at the edges. It is interesting to observe 
the interplay of this effect and the competing effect of 
the open boundary to lock in D^ ordering when is 
even (as demonstrated in Fig. [6]). This aspect of the 
VBS ordering is also important in view of DMRG stud- 
ies, where modified boundaries are often used.— "^"^ Here 
the boundary modification will simply be accomplished 
by excluding from the Hamiltonian the Q2 or Q3 terms 
with vertical bonds closest to an edge on every second 
row, as illustrated in Fig.[T4]in the case of Q2 terms. Re- 
sults obtained with only one of the edges modified will 
be compared with the case of both edges modified in the 
same way. 

The local variations of the VBS vector order parameter 
{Dx^ Dy) of the J-Q2 model were previously investigated 



14 



for 1/ X L lattices with all open edges<^ The formation of 
a vortex-like structure in the order parameter was noted. 
In the cases studied here, there is still translational sym- 
metry with period two along the y-axis and, therefore, a 
ID description of the order parameter as a function of the 
X coordinate suffices. The local x and y order parameters 
arc defined using the dimer operator in Eq. (|lll) : 



D,{x)= [{B.,{x,y))~^{B.,{x~\,y)) 

-i(i?,(x + l,y))](-ir, 
Dy{x) = [{By{x,y)) - {B,{x,y + 



(20) 



(21) 



These quantities are independent of the y coordinate (and 
an average is taken in the simulations to improve the 
statistics). A VBS angle 9{x) can also be defined. 



e{x) 



at an 



Dy{x) +Dy{x + 1) 
2DJx) 



(22) 



such that 9 — and 6* = tt for a fully x or y oriented VBS 
order, respectively. The reason for using the sum Dy{x) + 
Dy(x -I- 1) in the numerator under atan() is that an x- 
oriented column of bonds labeled by x is located between 
the j/-oriented columns at x and x + 1 (although such a 
detail of the definition of the local angle is not strictly 
important, and there are other equally good definitions 
giving the same result for large systems). 

Fig. [15] shows results for the Q2 model with only one 
modified edge. Oscillations in the bare dimcr expecta- 
tion value {Bx{x)) are present (top panel as in the case 
of the uniform edge in Fig. [6l In this case, however, the 
function is not reflection symmetric, due to the unequal 
left and right edges of the cylinder. The order param- 
eter Drc{x) is the largest at the edges. Away from the 
edge it decays toward a value at the center of the system 
which is somewhat smaller than the lockcd-in order pa- 
rameter previously extracted based on the data in Fig. [9] 
(which can be seen by analyzing data for several system 
sizes, not shown here). This is because the modified edge 
also leads to some amount of y order (middle panel of 
Fig. [151 and although this induced order decays rapidly 
when moving away from the modified edge it does not go 
away completely, even close to the opposite edge. 

The rather smooth decay of the y order to almost zero 
at the opposite edge can be explained as due to the open 
edge strongly favoring x ordering in its vicinity, even with 
the modification that breaks the y translational symme- 
try. The modified edge therefore induces both x and y 
order, i.e., the VBS angle (gH) is < 6* < tt/2. Since the 
second edge does not break the y translational symmetry 
explicitly, the x ordering can completely dominate there, 
leading to a very small (Dy) . The smooth transition from 
mixed x and y to almost pure x order is seen clearly in 
the VBS angle graphed in the bottom panel of Fig. [T5] 
Away from the edges, the total order parameter for large 
systems, D = {{Dx)"^ + (D^)^)^/^, approaches the value 
extracted for this model in the previous section. A max- 
imum in the angle develops with increasing size close to 
the modified edge. 
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FIG. 15: (Color online) Location dependent expectation value 
of the VBS order parameter of the Q2 model on 2L x L cylin- 
ders with the left (a; = 0) edge modified by the symmetry- 
breaking perturbation (inducing j/-oriented order) illustrated 
in Fig. 1141 The right edge is kept uniform. The top panel 
shows both the bare dimer expectation value {B{x)) and the 
dimer order parameter {Dx{x)) extracted from it according 
to Eq. (|20|l for L = 8. The middle panel shows the y order 
parameter defined according to Eq. (|21[l for L = 8, 16, and 
32. The bottom graph shows the VBS angle extracted from 
the X and y order parameters according to Eq. (f22|l . 



Fig. llGl shows results for systems with the y symmetry 
broken at both edges. Also in this case it would appear 
that both the x and y order parameters survive through- 
out the whole system in the thermodynamic limit. Con- 
vergence of both components as well as the angle at the 
center of the system is seen. The VBS angle here be- 
ing only slightly less than 7r/2 corresponds to an almost 
equal mixture of x and y order. 

In spite of the apparent convergence of the VBS angle 
to a value close to 7r/2 in Fig. [TBI the survival of both 
X and y order in the thermodynamic limit due to the 
modified edge is illusory. Since the VBS order is colum- 
nar, eventually, for very large systems, one would expect 
only X or only y order to survive. The explanation of 
the behavior seen is again the very large U(l)-Z4 cross- 
over scale in the Q2 model (as discussed in Appendix 
|A|. It is then interesting to look at the same quantities 
in the Q3 model, where there is no clear U(l) symme- 
try (as also shown in Appendix IX)) . i.e., the length-scale 
A is very short in this case. Results analogous to those 
in Fig. [THI for the Q2 model are shown in Fig. [T7| for 
the Qs model. In this case, one can see clearly how the 
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FIG. 16: (Color online) Location dependent expectation 
values of the x and y VBS order-parameter components, 
Eqs. PU)) and (PT|) . of the (32 model on 2L x L cylinders with 
both edges modified by a symmetry-breaking perturbation 
(favoring j/-oriented bond order). The corresponding VBS 
angle extracted using Eq. (|22p is shown in the bottom panel. 
The horizontal dashed line is at O = 7r/4 (corresponding to 
equal x and y order parameter parameters). 
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FIG. 17: (Color online) Same as Fig. [T5] for the Qs model. 
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FIG. 18: (Color online) The x and y components of the in- 
duced order parameter close to a modified edge of the Qs 
model on a 128 x 64 lattice. These data are the same as 
those shown in the top [x) and middle {y) panels of Fig. [T7] 
for smaller systems, but with the non-zero constant behav- 
ior at the center of the system subtracted off in the case of 
the X component. The lines are exponential fits, giving decay 
lengths 1.9 and 6.5 for the x and y components, respectively. 



y component vanishes with increasing system size away 
from the edges, while the x order stabilizes to a constant 
value. Since the x component is the surviving one, its 
approach to its bulk value should be governed by the 
standard VBS correlation length 5. The decay of the y- 
component should reflect A, however (since the presence 
of y order is due to the angular twisting of the order pa- 
rameter). This is a direct physical method to access the 
U(l) length-scale, providing an attractive alternative to 
studying the order-parameter distributions discussed in 
Appendix |^ 

The decays of the two componentss are analyzed quan- 
titatively for a larger system in Fig. [THl Excluding the 
points immediately adjacent to the edge, the decays are 
of almost pure exponential form (with an even-odd effect 
seen for the y component), giving ^ = 1.9 extracted from 
the X component and A = 6.5 from the y component. A 
similar analysis for the Q2 model (not shown here), based 
on systems with up to 128 x 64 sites, gives ^ « 25 (and 
A much larger still), but this estimate is not reliable be- 
cause the form of the decay is affected by the proximity 
to the critical point and is far from a pure exponential 
at the accessible distances. Larger system sizes are re- 
quired in this case, especially for extracting A, which is 
larger than 100 lattice constants according tho the analy- 
sis in Appendix [X] (perhaps being several hundred lattice 
constants). A systematic study of the divergence of the 
decay lengths of the J-Q3 model upon approaching the 
quantum-critical point will be presented elsewhere. 
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FIG. 19; (Color online) VBS order parameter and correlations 
lengths on infinite (Lx — >■ cxj) Q2 and Q3 cylinders of width 
Ly = 4, 6, 8, 10, 12. (a) The correlation length extracted using 
the dimer correlations with j/-oriented bonds as a function of 
Ly for those cylinders that have disordered ground states (the 
X correlation lengths are about 5 — 10% smaller), (b) The 
order parameter versus Ly, along with the corresponding 2D 
order parameters (shown with the horizontal lines). 



V. LONG CYLINDERS 

In the previous sections the 2D limit was approached in 
systems with fixed aspect ratio Lx/Ly. In principle the 
limit can also be accomplished with one of the lengths 
taken to infinity first, e.g., ^ 00 for fixed Ly and 
then Lj, — > 00. The behavior of the long-distance corre- 
lation functions, and, therefore, the squared VBS order 
parameter (D^), should not necessarily be expected to 
be smooth, however. Although VBS ordering amounts to 
breaking a discrete symmetry, and order can therefore, 
in principle, exist for any Ly in the infinitely long ID 
cylinder geometry, the survival of the order for small Ly 
is not guaranteed. Clearly, there will be enhanced fluc- 
tuations associated with the ID nature of these systems, 
which may destroy the ground-state order of a Hamilto- 
nian exhibiting long-range VBS order in the 2D limit. 

A well known system with discrete symmetry break- 
ing is useful for illustrating the potentially unsmooth ID 
to 2D cross-over: The Ising model with nearest-neighbor 
coupling Jz in a transverse magnetic field has a phase 
transition to an ordered (in the z spin direction) state 
at a critical value {h^/ Jz)c- On a ID linear chain the 
critical ratio is (hxj Jz)c = 1, while on the 2D square 
lattice it is (hx/Jz)c ~ 3.05i^ For an Lx x Ly lattice 
with Lx ^ 00 one can expect {hx/Jz)c to be a mono- 




FIG. 20; (Color online) VBS correlation functions, as de- 
fined in Eq. (|24p . for j/-oriented dimers in the Q3 model as a 
function of the separation in the a;-direction on cylinders in 
the Lx ^ 00 hmit. For Ly = 4, 6, fitted curves of the form 
C oc exp(— 2;/^)/a;" to the 2; > 4 data are also shown (with 
a ~ 0.5 in all cases). 



field 1 < hx/ Jz < 3.05, one can expect cylinders with 
small Ly to be disordered, while above some "critical" 
Ly the system will be ordered. One can expect the same 
kind of behavior of a 2D VBS as well, when restricting it 
to a finite cylinder, unless the 2D order parameter is ex- 
tremely large so that even the smallest cylinder remains 
in the ordered phase. 

In the discussion below, only cylinders of even Ly will 
be considered, so that the lattice is commensurate with 
columnar VBS order in both the x and y direction. J-Q 
models with odd Ly cannot be studied with the QMC 
method used here, because of sign problems arising due 
to geometric frustration of the spin interactions. 



A. Destruction of VBS order on cylinders 



tonic increasing function of Ly. Therefore, for a fixed 



As shown in Sec. IIIIBl the ground state of the pure 
model is strongly VBS ordered, the order parameter 
being at 70% of the maximum possible value. One might 
expect this to be sufficient for the order to be stable also 
on thin cylinders when Lx ^ 00. However, it turns out 
that such cylinders of width Ly = A and 6 are disor- 
dered, while for Ly ~ 8 and above the order parameter is 
already close to the 2D limiting value. For the pure Q2 
model, where the 2D order parameter is about 20% of the 
maximum value, no order was found on — > 00 cylin- 
ders with Ly up to 12. Larger widths were not studied 
due to prohibitively long computation times. The results 
for both models are summarized in Fig. [191 The results 
underlying these conclusions are discussed next. 

It is useful to define correlation functions that are av- 
eraged over the short (y) direction. The following func- 
tions, based on the definition PU)) of the elementary 
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dimer correlator, can be used to detect columnar VBS 
order with the bonds oriented either along the x or the 
y direction; 

Sd.,{x)^—Y,[Cd.{x,y) (23) 

-\Cdx{x~ l,y)- \Cdx{x + l,y)\, 

Sdy{x) ^ — Cay{x,y){-\)y . (24) 

Here it is appropriate to use periodic boundary condi- 
tions in both lattice directions. In order to achieve the 
limit ij. — > oo, aspect ratios L^jLy up to 32 were studied 
for Ly up to 12. 

In the Qa model, the y-dimer correlator Sdyix) ap- 
proaches a non-zero constant for large x when Ly > 8, 
as shown in Fig. [20l while for Ly = 4, 6 the correla- 
tions decays exponentially with distance. The behavior 
is not purely exponential but follow the form Sdy{x) oc 
a;~"exp(— cc/^), with a ~ 0.5. This form with a — 1/2 
is the Ornstein-Zernike (mean-field) form expected in a 
d = (1 + 1) dimensional system, where a = (d — l)/2. 
The correlation lengths extracted from fits to this form 
(with a regarded as a free parameter, to produce some- 
what better fits) are shown in Fig.fTWa'l. The x-oriented 
correlation function Sdx {x) is exponentially decaying for 
all Ly, i.e., these systems are purely y-ordered in the 
thermodynamic limit (as was also found in Scc. lIIIBl for 
periodic 2LxL systems when L — oo). For Ly = 4,6 the 
X correlation lengths are slightly smaller than the y ones. 
The y correlation lengths are graphed in Fig. [T9la). 

In Fig.[21]both the x and y correlation functions for the 
Q2 model are graphed for all even-width cylinders with 
Ly = 4, . . . , 12, along with fits to the exponential form 
discussed above. The y correlation length ^y is the larger 
one (about 5 — 10% larger than £^x) and is graphed versus 
Ly in Fig. [TWa). The correlation length grows roughly 
linearly with Ly for these cylinders. It would be interest- 
ing to go to even larger Ly to study the form in greater 
detail, and of course to find the threshold width for order- 
ing in this case (where presumably the correlation length 
should diverge, if one regards Ly as a continuous param- 
eter). The rather small correlation lengths for Ly up to 
12 suggest that it may be difficult to reach the critical 
width with QMC calculations at present. 

The destruction of the VBS order even on rather wide 
cylinders is surprising. In the Q3 model, judging by 
the decay of the x component of the order parameter in 
Fig. [THl the 2D correlation length is approximately 2 lat- 
tice constants. The lower width Ly = 8 for ordering on 
infinitely long cylinders is therefore roughly four times 
the correlation length. Moreover, related to the short 
correlation length, the 2D order parameter is as large as 
70% of the classical value. One might have expected such 
a system to be describable essentially in terms of classi- 
cal (orthogonal, hard-core) dimers with quantum fluctu- 
ations of the nature present in quantum dimer models. 
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FIG. 21: (Color online) VBS correlation functions, as defined 
in Eqs. (|23|l and p4|l . for x- (top panel) and y-oriented (bot- 
tom panel) dimers in the Q2 model as a function of the sep- 
aration in the x-direction on cylinders in the ^ 00 limit. 
Fits of the data for a: > 4 to the form S oc exp(— 2;/^)/a;" 
(woth Q ~ 0.5 in all cases) are shown as solid curves. 



It has been expected that a VBS under these conditions 
should be ordered even on narrow cylinders 1^ The re- 
sults obtained here suggest that the non-orthogonality of 
the singlets (the true quantum dimers) has a dramatic 
effects of reducing the order on cylinders, in contrast to 
this effect actually enhancing the dimmer-dimer correla- 
tions relative to those in corresponding dimer models in 
critical 2D systemsi^i^ On the other hand, to the au- 
thor's knowledge, quantum dimer models that order in 
the 2D limi1(2i have actually not been extensively stud- 
ied in long-cylinder geometry. Such studies would clearly 
be worthwhile, in light of the surprising results obtained 
here. 

In the Q2 model the correlation length should be in 
the range 20 ~ 30 (with, as already discussed above, the 
large uncertainty being due to the fact that system sizes 
-/j 3> C are needed to determine ^ accurately), and one 
can, thus, expect, roughly, Ly « 100 to be needed before 
ordering sets in on the cylinders in this case. 

One might speculate that the emergent U(l) symme- 
try could play some role in destroying the VBS order on 
the long cylinders. The local coarse-grained VBS order 
parameter (Dx,Dy) is an essentially isotropic 2D vector 
up to a large length scale A ^ ^1+° with a > (with the 
best estimate so far— being a = 0.2G±0.05). If the order 
parameter were truly a vector with isotropic angular fluc- 
tuations, long-range order on the ID Lx ^ 00 cylinders 
would be strictly prohibitedi^ The almost continuous or- 
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FIG. 22: (Color online) Boundary induced x and y compo- 
nents of the dimer parameter of the Qz model on 64 x 4 (top) 
and 128 x 8 (bottom) lattices. Both edges are modified to 
induce y order, as discussed in Sec. lIVI 



der parameter could then be argued to contribute to the 
loss of order. If so, one would expect a critical state to 
replace long-range order, however, of which there are no 
signs here — the VBS order decaying exponentially start- 
ing from short distances. There is no cross-over from 
a critical behavior, which might have been expected if 
almost U(l) symmetric angular VBS fluctuations were 
responsible for the destruction of long-range order. The 
role of emergent U(l) symmetry on cylinders is never- 
theless interesting and should be studied more systemat- 
ically in the future. 

Regardless of the exact relationship between the 2D 
correlation length and the ordering threshold on cylin- 
ders, the very short correlation lengths found in the Q2 
model (ranging from about ^ sa 2 for Lj, = 4 to ^ ss 8 
for Ly ~ 12) show the dangers of using the long-cylinder 
geometry for drawing conclusions about the presence or 
absence of VBS order in the 2D limit. Order likely ap- 
pears in the Q2 model, and probably in most models for 
which the existence of VBS order is under debate, for Ly 
far exceeding the maximum size that can currently be 
studied (especially if QMC methods cannot be used and 
DMRG would be the best choice of method). 

The above conclusions regarding ordered and disor- 
dered cylinders reached based on correlation functions in 
long periodic x Ly systems can also be confirmed by 
examining open-edge cylinders, in which a unique VBS 
can be locked in for even Lx (as discussed in the case of 
2Lx L cylinders in the preceding sections). Fig.l^shows 



results for longer Q3 cylinders in which the boundary 
perturbation inducing y order was also applied at both 
edges (as in Fig. [T3]). For Ly = 4, both order parameter 
components decay quickly away from the edges, while for 
Ly = 8 the y component stabilizes at the center of the 
system, at a value agreeing with that extracted on the 
basis of the correlation functions [shown in Fig. [TWa)]. 
Here, although the open edges favor x order more than 
the perturbations favor y order, the y component even- 
tually wins because that is the component favored just 
by having a finite Ly, and this effect scales with Lj-. In 
contrast, for the 2L x L cylinders with the same types of 
edges, it is the x component that survives in the thermo- 
dynamic limit, as seen in Fig. 1171 



VI. CONCLUSIONS AND DISCUSSION 
A. General summary and conclusions 

Several bench-mark results for the finite-size behavior 
of the VBS order parameter have been presented in this 
paper. The J-Q and pure Q models allowed investiga- 
tions of both strongly and weakly ordered ground states. 
The main general conclusion (which should be valid for 
VBS states in many systems) drawn from these studies is 
that even when the VBS order is relatively strong on the 
infinite 2D lattice (e.g., 10 — 20% of the maximum value 
attainable), results for small and moderate lattices (e.g., 
with up to hundreds of spins) can exhibit nearly critical 
behavior. The squared VBS order parameter then ap- 
pears to extrapolate to zero in the thermodynamic limit. 
In the J-Q model, this behavior can be traced to a rather 
large quantum-critical scaling regime around the critical 
value oi J/Q, where the behavior follows closely that ob- 
taining at a critical point. 

The extrapolation to infinite size may at first sight 
seem easier when symmetry-breaking boundaries are 
used (as is often done in the context of DMRG studies^) , 
so that the order parameter can be computed directly 
(having a considerably larger value than its square when 
the VBS order is not very strong). However, a small order 
parameter (10 — 20% of the maximum value in the VBS 
systems considered here) is very difficult to extrapolate 
accurately in this way, partially because the symmetry is 
not completely broken on lattices of size that can be stud- 
ied in practice. In particular, the emergent U(l) symme- 
try of the VBS order parameter implies that the com- 
ponent not locked by the boundaries can survive in the 
form of significant fluctuations up to very large system 
sizes, but this aspect of the ordering may be completely 
missed if one only examines the boundary-induced com- 
ponent of the order parameter. While this effect by itself 
would probably not lead to wrong conclusions regard- 
ing the presence or absence of VBS order, it is still im- 
portant for explaining results that would otherwise seem 
inconsistent with each other (e.g., when comparing the 
total squared order parameter and a direct boundary in- 
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duccd order parameter, as was done here in Sec. IIII Cp . 
The results presented here suggest that the best quantity 
for extrapolating the order parameter to infinite system 
size is the total (sum of the x and y components) long- 
distance correlation functions on L x L periodic lattices. 
Non-square lattices can lead to non-monotonic finite-size 
behavior. 

Some of the small-system behaviors pointed out here 
are generically well known and not limited to VBS order. 
There are also many examples of finite-size scaling of 
results for small lattices leading to wrong conclusions of 
the nature of the ground state. For example, in Refs. [9^ 
and [9^ a spin liquid ground state was claimed to exist 
in a 2D system of weakly coupled S = 1/2 Heisenberg 
chains. When QMC results for larger systems became 
available,— they showed a cross-over of the scaling and 
an asymptotic behavior in accord with a Neel state for 
any value of the inter-chain coupling. 

The additional complications due to emergent U(l) 
symmetryii^i^^ are more specific to VBS ordering. Open 
edges twist the vector order parameter {Dx,Dy) in ways 
which depends on the model and the nature of the edge. 
For a VBS there is no "neutral" edge; any boundary af- 
fects the ordering pattern in its neighborhood. While in 
the bulk VBS, in the thermodynamic limit, only one of 
the components can survive in a columnar state, at edges 
they can both be present. Due to the large length-scale 
of the cross-over from the U(l) symmetric order param- 
eter, both components can also survive in the interior of 
large systems. It would be interesting to study this phe- 
nomenon also in systems with a more complicated (larger 
unit cell) VBS order parameter. 

It should be noted that, although the concept of emer- 
gent U(l) symmetry of VBSs was developed in the con- 
text of deconfined quantum-critical points and has been 
confirmed in the case of J-Q models^i'^'^ this aspect of 
VBS order is most likely very general and manifested also 
in systems that arc not very close to such critical points 
(in some extended parameter space) — in 2D systems in 
which "angular" VBS fiuctuations are possible once the 
correlation length is several lattice constants or larger. 
The U(l) related boundary effects should be absent in 
cases where the angular fiuctuations are absent, e.g., in 
the case of staggered VBS statesj ^^'^°° 

For the purpose of detecting VBS order, an important 
aspect of the critical scaling is that, once a critical point 
has been identified, upward deviations from the power- 
law behavior, as seen in Fig. [T51 at J/Q2 ^ and 0.03 in 
the J—Q2 model, demonstrate an ordered state although 
this may not be apparent when carrying out extrapola- 
tions of the order parameter in 1/L (as in Fig. [T]). In 
general, in a model with some tunable parameter that 
can bring it into or out of a VBS state, it may be eas- 
ier to detect a phase transition than to extract the exact 
value of the order parameter close to such a point. On the 
one hand, many frustrated systems may have VBS states 
that are always only weakly ordered and, hence, close 
to a quantum critical point (or weakly first-order transi- 



tion) in some extended parameter space. Such systems 
should exhibit near-critical scaling on small lattices. On 
the other hand, if no critical scaling can be detected, and 
instead the order parameter correlation function decays 
exponentially fast with distance (or shows a tendency to 
decay faster than a power law) , one can rather safely con- 
clude that there is no VBS long-range order. Also with 
this approach, one can of course not expect to draw re- 
liable conclusions unless the system sizes are sufficiently 
large (and how large that is depends on the model). 

A striking behavior that may be particularly promi- 
nent in the case of VBS order was found here for lat- 
tices in the form of long cylinders, of size x Ly with 
Lx — > 00 and finite even Ly . In this geometry the order is 
unstable, and the system exhibits only short-range VBS 
correlations, until Ly exceeds some threshold that can be 
very large (perhaps 3-4 times the VBS correlation length, 
according to results for the Qt, model). Long cylinders 
arc therefore not ideally suited for determining the na- 
ture of the 2D state in which VBS order is a possibility 
(systems with a small fixed L^/Ly normally being bet- 
ter). In particular, the method of positively confirming 
a Z2 spin liquid by the absence of of order on even-Ly 
systems is not applicable in the "yes-no" sense proposed 
in Ref. [H. Instead, the finite-size behavior has to be 
tracked as in any other extrapolation method. The cor- 
relation length as a function of even Ly should converge 
for a spin liquid and diverge for a VBS, as in Fig. [T9l 
but it may not be easy in practice to determine which of 
these behaviors applies. 



B. Comment on the possibility of a spin-liquid 
state in the J1-J2 Heisenberg model 

One motivation for the present study was to pro- 
vide guidance on detecting VBS order — or, alternatively, 
showing the absence of such order — in calculations for 
frustrated 2D models. The lattice sizes reachable for such 
systems with unbiased calculations, primarily using the 
DMRG methodr'^ are stiU very limited. Methods based 
on tensor-product states^^i^i^ beyond matrix-product 
states (which are closely related to the DMRG scheme), 
are still typically too much affected by various truncation 
errors and approximations to be considered completely 
unbiased. The following discussion will therefore be pri- 
marily aimed at DMRG calculations, although many of 
the issues would apply more generally. 

The issues raised here have particular relevance in the 
context of a recent DMRG study of the J1-J2 Heisenberg 
model on the square lattice^i^ Several different ways of 
analyzing VBS correlations were argued to consistently 
show the absence of VBS order and positively confirm 
the properties of a Z2 spin liquid. However, many of the 
results presented can also be explained by a VBS state, at 
least in some part of the non-magnetic phase, according 
to the results obtained here. The key points supporting 
this view are summarized next. 
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In Fig. 3 of Ref. @, second-order polynomial fits to the 
VBS order parameter for 2LxL cylinders with L < 10 are 
shown. The fact that these fits extrapolate to negative 
values in the thermodynamic limit was taken as evidence 
for the absence of VBS order. However, this kind of be- 
havior is also observed for the Q2 model on small lattices, 
as seen in Figs. [7] and [S] of the present paper, even though 
the order parameter of this model is as large as 20% of 
the maximum possible value. If VBS order exists also in 
the non-magnetic phase of the J1-J2 Heisenberg model, 
one should not expect it to be very strong. Therefore, 
the finite-size behavior seen in Fig. 3 of Ref. is at least 
qualitatively what would be expected even if the state is 
a VBS. It should be noted that the fact that the fitted 
functions extrapolate to negative values is in itself a clear 
sign of the chosen functional forms not being correct, as 
the squared order parameter cannot be negative. Thus, 
there must necessarily be a cross-over to a different form 
for larger systems — either to a pure 1 /Li^ form, if there is 
no long-range order, or to an exponentially rapidly con- 
vergent form tending to a non-zero value. The results for 
small systems cannot distinguish between these different 
asymptotics. 

The finite-size extrapolation issues may clearly also af- 
fect the determination of the transition point between 
the Neel antiferromagnet and the non-magnetic state at 
g = J2/J1 ~ 0.4 (while the transition point into the 
stripe antiferromagnet at g = J2/ Ji « 0.6 is much easier 
to extract due to it being clearly first order). The tran- 
sition point g « 0.41 was determined in Ref. @ based on 
extrapolations of the Neel order parameter (M^) using 
second-order polynomials, and these should be affected 
by similar problems as those pointed out here for the 
VBS scaling (and it is also well known that polynomials 
higher than second order have to be used to extrapolate 
Neel order correctly based on small systems, even in the 
strongly order Heisenberg model^SiSi). The Neel order 
should therefore survive up to somewhat larger g values. 
Thus, at g = 1/2, on which most of the analysis of the 
VBS scaling was focused in Ref. 0, the system may be 
rather close to the transition point. If VBS order exists 
in the nonmagnetic phase, it would therefore likely be 
very weak at this point. In Fig. 3(a) of Ref. 0, the max- 
imal value of the order parameter (Dy), at g just below 
0.6, is close to the values for the Q2 model in Fig. |9]of 
the present paper. Thus, if the J1-J2 model has VBS 
order, its peak value should be about 10 — 20% of that of 
a perfect columnar state. It would be better to analyze 
the VBS correlations closer to the maximal value, where 
the extrapolation problems are minimized. 

As discussed in Sec. IIIIDl in systems where there is a 
quantum phase transition into the state of interest, the 
best way to deduce the nature of that state may be to 
first carefully examine the phase transition. If there is 
critical scaling, deviations from the power-law form of 
the order parameter away from the critical point can be 
a good signal of long-range order. However, as seen in the 
scaling plot for the near-critical J-Q2 model in Fig. [T^l 
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FIG. 23: (Color online) Finite-size scaling of the squared VBS 
order parameter (a) and staggered magnetization (b) calcu- 
lated on the central L x L square of cylinders of size 2L x L. 
DMRG Results for the J1-J2 model at g — 0.50 and 0.56, from 
Figs. 2(a) and 3(a) of Ref. 12, are compared with QMC results 
for the J-Q2 model at its critical point, [J/Q2)c ~ 0.0447, and 
at J = 0. In (a) the VBS y component of the J1-J2 model 
and the x component (the larger component) of the J-Q2 are 
shown. The line drawn close to the J/Q2 ~ 0.0447 points has 
slope —1.27, corresponding to the critical exponent r; = 0.27 
(as in Fig. I13p . and that going through the g = 0.50 points 
has slope —1.8. The dashed line has slope —2, corresponding 
to the expected asymptotic behavior in a non-VBS state. In 
(b) both lines have slope -1.27 (77 = 0.27). 

if the accessible system sizes are only up to L w 10, 
even a system in which the VBS order parameter is as 
large as 20% of the maximum value may in practice not 
be distinguishable from a critical system when analyzing 
the order parameter fiuctuations. If the non-magnetic 
state of the J1-J2 Heisenberg model also has long range 
order, then one should expect a similar behavior. 

Re-plotting the g = 0.5 and 0.56 data for the VBS 
y component of Fig. 3(a) of Ref. on a log-log scale, 
one can indeed observe behaviors close to power laws, as 
shown in Fig. [53Ja). In the same graph data for the J- 
Q2 model at J = and Jc ~ 0.0447 are also graphed. 
In this case the x component of the order parameter is 
shown, which, as seen in Fig. [51 in this system is larger 
than the y component and is the one surviving in the 
thermodynamic limit. In the J1-J2 model it is instead 
the X component that is somewhat larger i^SL 

The comparison of the two models is complicated by 
the fact that the average induced x order was subtracted 
in the definition used in Ref. 0. That induced order is 
very small, however fiSi unlike what it is in the Q2 model 
(which, may indicate that the VBS order, if it exists in 
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the J1-J2 model, is y-oricnted on the cyhndrical 2L x L 
systems, as was also noted in Ref. @). 

For the open-edge 2L x L cylinders used in Fig. E^ a) , 
the J-Q2 results do not exhibit quite as good scaling as 
in the case of the periodic L x L systems in Fig. [131 
but for large systems the behavior is still consistent with 
an exponent 77 « 0.3. The J1-J2 results for g = 0.5 
follow a different behavior, however, decaying as L^" 
with a ~ 1.8. This is quite close to a = 2, which is 
expected deep inside a non-VBS phase. For g = 0.56 
the data for the larger sizes deviate significantly upward 
from the g ~ 0.5 points and cannot be fitted very well to 
a power law. Tthe slope on the log- log scale is ~ —1.52 
for a line drawn through the L = 8 and 10 points, but the 
data for smaller systems fall above the fitted line, showing 
a flattening out with increasing size. The reduction of 
the rate of decay is opposite to the expectation for a spin 
liquid and an indication that the system is VBS ordered 
in the infinite-size limit. 

The behavior at g = 0.5 is puzzling. Since the VBS 
order parameter here follows quite close to the form ex- 
pected in a spin liquid, one may conclude that this is 
what it is, and the deviations from the '-^ form arc 
due to remaining size effects (i.e., the system size is not 
yet much larger than the correlation length). A possibil- 
ity suggested by the behavior observed in Fig. [23] is that 
the J1-J2 model has a spin liquid phase following the 
Neel phase above g « 0.4, followed in turn by a VBS at 
larger g (since the g = 0.56 results seem more indicative 
of weak VBS order). Another possibility is that there is 
no spin liquid, but the Neel- VBS transition takes place 
at g significantly larger than previously believed, so that 
g ~ 0.5 would actually still be inside the Neel phase. 
Looking at the raw data for the sublattice magnetiza- 
tion in Fig. 2(a) of Ref. 0, it appears that this possibility 
cannot be ruled out (considering again also the fact that 
the second-order polynomial fits used should lead to an 
under-estimation of the critical g where the Neel order 
vanishes). The behavior of the triplet gap in Fig. 2(b) 
seems to go against this scenario, however, although the 
way the gap was extracted, by targeting higher states ob- 
tained while keeping the edges in the ground state, may 
lead to strong corrections to the gap scaling. 

To investigate possible near-criticality in the Neel order 
parameter, the results from Fig. 2(a) of Ref.^for (M^) at 
g = 0.5 are re-plotted on a log-log scale in Fig. l^^ b). In- 
terestingly, the behavior follows closely a power law, with 
an exponent 77 very similar to that of the J-Q2 model. 
This could indicate that the transition out of the Neel 
state indeed takes place close to g = 0.5 and is in the 
same universality class as the J-Q model. Note that, 
within the deconfined quantum criticality theory,— this 
kind of criticality of the magnetic order would not neces- 
sarily require that the VBS order emerges at this point 
as well, because the exponents associated with the Neel 
order parameter are not affected by the VBS (since the 
operator causing the VBS order is dangerously invari- 
ant). Clearly there is not sufficient data here to make 
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FIG. 24: (Color online) VBS y order parameter component 
induced by edges modified to break the y translational sym- 
metry (as explained in Fig. I14p in the Q2 model. Here r is 
defined as the distance from the second column of spins away 
from the edge, since the edge modification extends to this lo- 
cation. The lines show exponential fits, with decay lengths 
1.8 {Ly = 4), 3.2 {Ly = 6), 4.8 {Ly = 8), and 6.6 {Ly = 10). 

any firm conclusions about this scenario of a Neel to spin 
liquid transition, possibly followed by a subsequent liq- 
uid to VBS transition at higher g, but the behavior is 
intriguing and deserves further tests. 

An important aspect of the analysis of Ref. 0, cited 
as positive evidence for a Z2 spin liquid, is the behav- 
ior of the order parameter on infinitely long cylinders. 
There is an even-odd effect that had previously been 
found in liquid states of quantum dimer models)^ For 
odd Ly and even 00, an x-oriented order pa- 

rameter cx exp(— Ly/^j,) is induced because of geometric 
frustration effects. For even Ly no order is observed at 
all, regardless of the type of VBS (horizontal or vertical 
columns) favored by the edges. Unfortunately, odd-Lj, 
J-Q cylinders cannot be studied with the QMC method 
used here, because of sign problems. However, based on 
the results presented here for even Ly it is already clear 
that this kind of test for a Z2 spin liquid may not be that 
useful in practice, because VBS order does not exist on 
the infinitely long cylinders (for Ly up to some critical 
width that can be expected to be inaccessible in practice 
for systems that are weakly to moderately ordered in the 
2D limit). In Ref. @ it was implicitly assumed that any 
system with 2D VBS order will exhibit such order also 
on long thin cylinders. 

It is also interesting to note that the induced order 
parameter as a function of the distance from a modi- 
fied edge of systems in the L^ ~> 00 limit is very similar 
in the J1-J2 and J-Q2 models. Fig. [24| shows results 
for the pure Q2 model on cylinders of width 4 — 10 in 
which the edge has been modified to break the y transla- 
tional symmetry, as discussed in Sec. IIVI and illustrated 
in Fig. [T4| Here cylinders with aspect ratio L^/Ly = 16 
were used (which is large enough to accurately represent 
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the Lx ^ oo limit). The edge- induced x and y order 
parameters both decay exponentially, with very similar 
decay lengths that are also close to the correlation lengths 
graphed in Fig. [TO] (obtained from correlation functions 
on systems with all periodic boundaries). The y decay 
lengths are always marginally larger. For i = 6 and 8, 
the decay lengths are about 1.5 times those in the J1-J2 
model at g = 0.5, for which data were shown in Fig S6(b) 
of Ref . 1 

As discussed above and seen clearly in Fig. [23l the 
VBS order parameter is likely significantly suppressed 
at g = 0.5 relative to what it is close to its maximum 
in this model (which appears to be a bit above 0.56). 
One can therefore expect to see decay lengths as large 
as those in the Q2 model for larger g (close to 0.6). The 
rapid decay was in Ref. H interpreted as the system being 
insusceptible to VBS ordering even in the presence of, 
at first sight, very favorable conditions for inducing it. 
Again, when analyzed in light of the known physics of 
the J-Q2 model, the results cannot be distinguished from 
those of a rather substantially ordered VBS. It would be 
illuminating to have J1-J2 data for Ly > 8, to see if the 
decay length continues to grow or saturates. 

In Ref. 1^ the size dependence of the entanglement en- 
tropy was also used as positive evidence of a Z2 spin liq- 
uid. It would be very interesting to compute this quantity 
also for the J-Q models. It is clear that the non-trivial 
aspects of the VBS fluctuations could lead to behaviors 
not predicted in the strong- VBS limit. Since the system 
on small lattices and cylinders resembles a spin liquid, 
it would not be surprising if the corrections to the area 
law of the entanglement entropy are also similar, up to 
some large size where the true asymptotic VBS behavior 
sets in. QMC calculations of the entanglement entropy 
of the J-Q models will be carried out in future studies, 
using the recent developments of methods to study the 
Renyi versions of the entropies This should clarify 

whether the constant deviation from the area law cited 
in Ref. is really unique to Z2 spin-liquids, or whether 
they can also appear (for lattices of practically reachable 
size) in weakly ordered VBS states. The scaling of the 
entanglement entropy at a deconfined quantum-critical 
point is also of interest heieJ^ 

The conclusion reached from the above comparisons of 
results for the J1-J2 model and the J-Q models is that 
they exhibit rather similar behaviors, and, therefore, a 
VBS ground state of the J1-J2 cannot be excluded. Some 
of the J1-J2 results may also be consistent with a Z2 spin 
liquid at g 0.5, but the point to note here is that most 
of the results presented so far do not favor that kind 
of state over a VBS state. In particular, the claimed 
positive signals for a Z2 spin liquid are also seen in the 
confirmed VBS state of the J-Q models. If anything, the 
very similar behaviors seen in the near-critical Q2 model 
and the J1-J2 models should tilt the balance further in 
favor of VBS order for g = J2/J1 close to 0.6. The be- 
havior at g = 0.5 is very intriguing and not consistent 
with a near-critical VBS of the same kind as in the J-Q 



models. It would be very useful to analyze the VBS and 
magnetic correlations further in this case, preferrably on 
larger lattices. 

It would also be good to know in greater detail the 
effects of truncation (the number of states kept) in the 
DMRG calculations. The error « 10^^ in Ref. prefers to 
the missing weight in the density matrix. One can expect 
the errors in the wave function to be approximately the 
square-root of this error^^iS but exactly how much the 
VBS correlations are affected, especially for the largest 
systems, is not entirely clear. 



C. Remarks on other potential spin liquids 

The results presented here also are relevant to stud- 
ies of the kagome Heisenberg model, for which DMRG 
studies also have indicated a spin liquid statei^i^ A VBS 
is another candidate state,— which is not easy to ex- 
clude if the ordering is weak (which should be expected, 
if this kind of order is present). Since the most likely 
VBS patterns in this case are much more complicated 
than the columnar state of the J-Q models discussed 
here (with the most likely candidate states having 12- 
or 36-spin unit cells) , it is not possible to relate results in 
the same close manner as done above in the case of the 
J1-J2 model. Nevertheless, the issues pointed out here 
should be considered also when analysing the kagome 
system, in particular on long cylinders. It would be very 
desirable to reach larger x Ly lattices with the as- 
pect ratio Lx/Ly kept fixed, although this seems difficult 
at present. It would also be good to push calculations 
based on the multi-scale entanglement renormalization 
ansatz (MERA)2i to higher precision. Such a calculation 
had previously seemed to confirm the VBS with 36-site 
cell proposed earlier based on other techniques,— but 
the energy reached was not as low as that found with 
DMRG^i^ and exact diagonalization.i^^ 

The analysis and arguments presented in this paper 
also suggest that it would be very useful to add to the 
nearest-neighbor Heisenberg exchange some term that 
favors one of the VBS states proposed previously, and 
to study the phase transition out of this ordered state. 
Longer-range couplings may work, but some interaction 
similar to the multi-spin Q terms discussed here could be 
even better suited for inducing the desired type of VBS. 

Spin liquid states have recently also been claimed to 
exist in electronic Hubbard models and frustrated spin 
models on the honeycomb lattice. ^^^^-^SI For the Hub- 
bard model, 2D lattices with up to hundreds of sites were 
usedJ^ The VBS correlations in this case decay very 
rapidly with distance, and the system does not seem to 
exhibit the kind of problematic scaling issues pointed out 
in this paper. On the other hand, work on effective spin 
models constructed to capture the putative spin liquid 
state have not so far been conclusive . ^'^'^"^" ■^ Also here 
it would be useful to extend the models in such a way 
that a VBS phase transition can be studied. The VBS 
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FIG. 25: (Color online) VBS order parameter distribution 
P{Dx, Dy) in the Q3 model on periodic L x L lattices with 
L = 12 (left) and L = 24 (right). The size of both squares 
corresponds to the full space of possible values of the com- 
ponents Dx,Dy £ [--Dmax, Anax], where Anax = 3/8 (for a 
perfect columnar VBS). 
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FIG. 26: (Color online) Size dependence of the columnar 
anisotropy weight, defined in Eq. (|A1|) . of the VBS order pa- 
rameter distribution in the Qs model. 



should then be the one to which the "bare" honeycomb 
model is the most susceptible (which may in itself not be 
easy to determine in this case). 



D. Bench-mark challenge 

Finally, as a challenge to DMRG, tensor-product, and 
MERA techniques, it would be very interesting and use- 
ful to see these methods applied to J-Q models as well. 
Comparing with the known phase diagram and critical 
behavior extracted on the basis of unbiased QMC simu- 
lations would be a very good test of the capabilities of 
these methods to capture non-trivial ground states and 
quantum phase transitions. If the outcome is positive, 
it may be very useful to systematically investigate the 
behavior when frustration is added to this model, as was 
recently done in an exact diagonalization study of a 2D 
model combining the Q2 interaction with the frustrated 
Ji- J2 Heisenberg modelJ^ 
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Appendix A: U(l)— Z4 cross-over of the VBS 
symmetry in periodic systems 

The emergent U(l) symmetry of a columnar VBS in 
the neighborhood of a critical point can be characterized 
by the probability distribution P{Dx, Dy) generated in 
QMC simulations on periodic L x L lattices. A system- 
atic study aimed at extracting the scaling of the U(l)- 
Z4 cross-over length A was presented in Ref. |4l|. Here 
additional results for the pure Q2 and Q3 models will 
be presented in order to facilitate comparisons with the 
boundary effects discussed in the main text. Specifically, 
it will be shown that the lack of D^-Dy symmetry on 
2Lx L lattices, as seen in Fig. 0] for the Q3 model for all 
system sizes, is matched by a clear Z4 symmetric order 
parameter on all L x L lattices. Conversely, the sym- 
metry seen for the Q2 model on large lattices in Fig. |3] 
is consistent with only very small deviations (barely de- 
tectable) from U(l) symmetry on L x L lattices with L 
as large as 128. 

In the projector QMC simulations, each generated con- 
figuration is associated with a pair of order parameters 
{Dx,Dy), which are matrix elements of the correspond- 
ing operators defined in Eqs. and computed 
in the valence bond basis. These matrix elements are 
of the form 37i/4A^, where n is an integer in the range 
[— A^/2, A^/2], with the extremal values corresponding to 
both the bra and ket state (making up the transition 
graph) having the same perfect columnar pattern of va- 
lence bonds of length one lattice constant. The his- 
togram P{Dx, Dy) is constructed based on these matrix 
elements. 

Fig. [55] shows results for the Q3 model for L = 12 and 
L = 24. In this model the histogram P{Dx, Dy) exhibits 
a distinct four-fold symmetry even for the smallest sys- 
tems (also smaller than L = 12, not shown here, where 
the discreteness of the distribution function also becomes 
apparent). The four peaks sharpen with increasing lat- 
tice size, and above some size the suppression of the 
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FIG. 27: (Color online) VBS order parameter distribution 
P{Dx, Dy) in the Q2 model on periodic L x L lattices with 
L = 64 (left) and L = 128 (right). The size of both squares 
corresponds to 10% of the maximum value Dmax/lO of the 
components, D^^^Dy £ [-Dmax, -Dmax], where Anax = 3/8 
(for a perfect columnar VBS). 

weight between the peaks severely impedes QMC fluctu- 
ations between the peaks. In fig. [551 the visibly different 
weight in the four peaks (with the right peak having the 
smallest weight) is a consequence of this rarity of "in- 
stanton" events between the peaks (i.e., the simulations 
"get stuck" in one quarter of the configuration space). It 
should be noted that this very slow simulation dynamics 
of the VBS order parameter does not affect the estimate 
of the total squared order parameter {D"^) and most other 
physical quantities of interest. 

The degree of Z4 symmetry of the order parameter can 
be quantified by the function 

Dy 

where 4>xy is the angle corresponding to the point 
{Dx,Dy). While this function (and the underlying prob- 
ability distribution) is not a physical observable, in the 
sense that it is not a bona fide quantum mechanical ex- 
pectation value, it nevertheless reflects the fluctuations of 
the VBS order parameter and can be used to characterize 
the the U(l)-Z4 cross-over. 

Results as a function of L for the Q3 model are shown 
in Fig. [551 Here the convergence W4 — >■ 1 when L — > 00 
is apparent, as would be expected for a columnar VBS 
in the thermodynamic limit. In principle the curve 
Wi{L) could be used to define the length A, e.g., us- 
ing W4(A) = 1/2, but there is clearly an arbitrariness in 
choosing the particular number. For studying the scaling 
of A when some parameter of the Hamiltonian is changed 
(e.g., J/Q3) this ambiguity does not matter. In Ref. l4ll 
curves W4{L) for different coupling rations were analyzed 
using standard finite-size scaling techniques, with the re- 
sults that A grows slightly faster than the correlation 
length; A - with a « 0.2. 

Comparing with the behavior of the squared order pa- 
rameters in Fig. [31 it can be noted that (-D^) approaches 
(and (Dy) tends to a non-zero value) very quickly above 



L 20, which is approximately where W4{L) = 1/2 in 
Fig.[26l On the other hand, the decay of the edge-induced 
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FIG. 28: (Color online) Angular distribution of the VBS order 
parameter of the Q2 model for system sizes L — 32, 64, and 
128. To improve the statistics, these results were obtained 
by symmetrizing the distributions using the expected 90° ro- 
tational symmetry. The jaggedness of the curves (especially 
for L = 32) is due to the discreteness of the allowed {D^, Dy) 
values (with A'' possible values for each component). 



y component of the order parameter in Figs. [T7l and [18] 
(where the system far from the edge has only x order) 
gives a length « 6.5, which could also be taken as a prac- 
tical definition of A. This length corresponds to W4 w 0.1 
in Fig. [H 

In contrast to the Q3 model, in the Q2 model no clear 
Z4 symmetry is visible in P(Dx, Dy) up to systems as 
large as L = 64 and 128, as shown in Fig. [571 These 
histograms are ring-shaped, although for L = 128 the 
weight is not evenly distributed because of lack of suf- 
ficient QMC statistics. The VBS angle fiuctuates very 
slowly in simulations of large systems and very long runs 
are required in order to obtain symmetric distributions. 
The data shown arc based on w 3.5 x 10® Monte Carlo 
sweeps for L ~ 64 and 8 x 10^ for L ~ 128 (which re- 
quired more than 10** CPU hours in both cases). By 
symmetrizing the distributions using 90° rotations, one 
can still detect small deviations from perfect U(l) sym- 
metry, as shown in Fig. [5H1 The peak positions again 
correspond to a columnar state. 

Note that in Fig. [571 the ring for L = 128 is consid- 
erably thinner than for L ~ 64, with the radius (the 
location of the maximum or average weight) remaining 
almost unchanged. This reflects an expected reduction 
of the fluctuations of the magnitude of the VBS order 
parameter with increasing system size. 

Based on these results, the cross-over length-scale A 
for the Q2 model should be ^ 128, which explains why 
both order-parameter components are essentially equal 
for the largest systems in Fig. [51 
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